purse.models.new
LiaC
2018-06-03
Model selection for Ch 1 - eelgrass and sand flat data
- Decided we should include site as a random effect for all purse seine models because the # of sites is well above 5 and so we deemed it important to acknowledge the nested nature of our sampling design.
- Early on in model selection we determined that our data has a negative binomial distribution. All models started with Poisson distributions and all were overdispersed, which was corrected by the neg binomial distribution. All models here use negative binomial distribution with a log link.
- We found that the relationship between abundance and Julian day is commonly better described as quadratic (vs. linear), and have incorporated this into the global model (i.e including Julian day as two parameters: Jday + Jday^2). In some cases this term was removed from the global model if it did not converge (see purse.chinook.Rmd)
- Habitat variables were tested for collinearity, and through this process we narrowed them down to Leaf Area Index and Mean Turbidity only to include in the global model. These two variables still do covary, but make biological sense as separate variables.
4.5) Site E5 (South Ferry) was missing measurements for leaf area index, so the values from site E6 (Point Roberts) were used. Site E7 (near coal port in intercauseway) was missing values for turbidity, so the mean of sites E3 and E4 (Intercauseway N and S) were used. Archipelago values for leaf area index were used for sites E3, E4, and E7. All eelgrass measurements were set to 0 for sand flat sites. Mean turbidity values were calculated for all sites (including sand flat)
- We used AICc to rank all models in dredge, as it penalizes the strength of likelihood more for very small sample sizes, and is more adaptive to moderate sample sizes, so that we could apply the same information criterion across all models (Harrison et al 2018; Brewer et al 2016) – could explain better…
- We used model averaging following a delta AICc of less than 4, [[and followed the nesting rule to remove any artificially added models that were represented by smaller models with better AICc scores (Harrison et al 2018).]]- didnt do this because weights were low enough anyway so arbitrary. Kept at delta 4 because weights dropped off considerably before even that point. The choice of cuttoff is somewhat arbitrary, but within the accepted standard of ranges to incorporate enough uncertainty while also having some selectivity. From Jenn:
Hi Lia,
I don’t have a reference that explicitly says that you should use a delta AIC cut off of 4, but my decision to use 4 as a cutoff was mainly based on two things:
- In Burnham and Anderson’s 2002 book they say “… an alternative rule of thumb for an approximate 95% confidence set on the K-L best model is the subset of all models having delta less than or equal to some value that is roughly 4-7” (p. 170), which seems to suggest that any value within that range would be a reasonable choice. I also found that I had a lot more models in the top model set when using a higher delta AIC cut off value (e.g. 7), so 4 seemed like a better choice.
- A lot of the stats for my paper are modeled after the stats used in another coral reef-related paper by Emily Darling (https://link.springer.com/article/10.1007/s00338-017-1539-z), where they also used a delta AIC cut off of 4 (although they don’t provide a reference for it).
Chinook purse seine model
response = average abundance of Chinook salmon [per purse seine site]
This is 1st of 4 separate purse models for (1)Chinook, (2)chin, (3)other migratory species and (4)resident species
Load all data - standardize variables, create subgroups from purse seine data
# load Purse seine data aggregated by species/set; contains unidentified
# larval
purse <- read.csv("/Users/Lia/Documents/Git/Fraser-salmon/all.data/purse.catch.csv")
# get summary stats for each model parameter for appendix
summary(purse)## Year Date Temp.surf DO.surf
## Min. :0.0000 6/2/2017 : 139 Min. : 7.49 Min. : 43.1
## 1st Qu.:0.0000 5/19/2017: 111 1st Qu.:12.51 1st Qu.:101.0
## Median :0.0000 7/15/2017: 99 Median :14.06 Median :112.2
## Mean :0.4481 6/3/2017 : 90 Mean :13.90 Mean :110.1
## 3rd Qu.:1.0000 6/20/2017: 87 3rd Qu.:15.31 3rd Qu.:119.1
## Max. :1.0000 24-Jun-16: 81 Max. :20.41 Max. :156.5
## (Other) :1620
## DOmg.surf pH.surf Sal.surf Temp.bot
## Min. : 4.04 Min. :7.780 Min. : 0.55 Min. : 6.19
## 1st Qu.: 9.25 1st Qu.:8.060 1st Qu.:12.96 1st Qu.:11.84
## Median :10.17 Median :8.250 Median :23.02 Median :13.22
## Mean :10.08 Mean :8.268 Mean :20.63 Mean :13.33
## 3rd Qu.:10.99 3rd Qu.:8.520 3rd Qu.:28.80 3rd Qu.:14.79
## Max. :14.56 Max. :9.000 Max. :33.92 Max. :20.22
##
## DO.bot DOmg.bot pH.bot Sal.bot
## Min. : 45.3 Min. : 4.13 Min. :7.02 Min. : 1.31
## 1st Qu.:100.6 1st Qu.: 9.16 1st Qu.:8.03 1st Qu.:21.11
## Median :110.0 Median : 9.79 Median :8.22 Median :27.17
## Mean :111.1 Mean :10.09 Mean :8.23 Mean :24.61
## 3rd Qu.:121.5 3rd Qu.:11.19 3rd Qu.:8.41 3rd Qu.:29.97
## Max. :177.3 Max. :16.85 Max. :9.20 Max. :34.05
##
## J.date SIN_TIME COS_TIME Round
## Min. : 79.0 Min. :-0.99771 Min. :-0.9999 Min. : 1.000
## 1st Qu.:133.0 1st Qu.: 0.01075 1st Qu.:-0.9765 1st Qu.: 5.000
## Median :154.0 Median : 0.47276 Median :-0.8644 Median : 6.000
## Mean :161.3 Mean : 0.32371 Mean :-0.6879 Mean : 6.173
## 3rd Qu.:182.0 3rd Qu.: 0.75370 3rd Qu.:-0.5185 3rd Qu.: 8.000
## Max. :286.0 Max. : 0.99999 Max. : 0.2102 Max. :10.000
##
## Habitat Site Set
## Eelgrass :1435 E5 :275 Min. :1.000
## Sand flat: 792 E3 :265 1st Qu.:1.000
## E2 :255 Median :2.000
## E4 :235 Mean :2.007
## E1 :205 3rd Qu.:3.000
## SF6 :178 Max. :4.000
## (Other):814 NA's :1
## Species Life_stage abundance
## Shiner surfperch :329 0 : 82 Min. : 0.00
## Three-spined stickleback:320 adult :884 1st Qu.: 1.00
## Starry flounder :293 juvenile:962 Median : 2.00
## Bay pipefish :147 NA's :299 Mean : 19.77
## Chinook :131 3rd Qu.: 6.00
## Surf smelt : 84 Max. :1698.00
## (Other) :923
## class round month
## 0 : 82 Min. : 1.000 April :357
## migratory: 370 1st Qu.: 5.000 July :341
## resident :1717 Median : 6.000 June :538
## NA's : 58 Mean : 6.078 March :153
## 3rd Qu.: 7.500 May :656
## Max. :10.000 October : 80
## September:102sapply(purse, sd)## Year Date Temp.surf DO.surf DOmg.surf pH.surf
## 0.4974146 14.0883494 2.5111848 15.2915703 1.4735729 0.2599412
## Sal.surf Temp.bot DO.bot DOmg.bot pH.bot Sal.bot
## 9.4074621 2.4150753 17.5070536 1.6823201 0.2550196 7.5839695
## J.date SIN_TIME COS_TIME Round Habitat Site
## 44.7545523 0.5426650 0.3573595 2.2216665 0.4788129 3.8004176
## Set Species Life_stage abundance class round
## NA 12.4404013 NA 80.1448604 NA 2.0653150
## month
## 1.6950018# eelgrass summary
eelgrass <- subset(purse, Habitat %in% "Eelgrass")
summary(eelgrass)## Year Date Temp.surf DO.surf
## Min. :0.0000 6/2/2017 :139 Min. : 7.49 Min. : 57.5
## 1st Qu.:0.0000 5/19/2017:111 1st Qu.:12.23 1st Qu.:105.0
## Median :0.0000 7/15/2017: 99 Median :13.85 Median :113.6
## Mean :0.4544 6/20/2017: 87 Mean :13.64 Mean :111.9
## 3rd Qu.:1.0000 24-Jun-16: 81 3rd Qu.:15.04 3rd Qu.:122.7
## Max. :1.0000 10-Jul-16: 77 Max. :18.85 Max. :156.5
## (Other) :841
## DOmg.surf pH.surf Sal.surf Temp.bot
## Min. : 5.530 Min. :7.780 Min. :10.38 Min. : 7.07
## 1st Qu.: 9.190 1st Qu.:8.110 1st Qu.:22.81 1st Qu.:11.63
## Median :10.160 Median :8.300 Median :27.36 Median :13.06
## Mean : 9.965 Mean :8.313 Mean :25.56 Mean :13.06
## 3rd Qu.:10.880 3rd Qu.:8.550 3rd Qu.:29.64 3rd Qu.:14.65
## Max. :14.190 Max. :9.000 Max. :33.92 Max. :17.53
##
## DO.bot DOmg.bot pH.bot Sal.bot
## Min. : 64.9 Min. : 5.96 Min. :7.02 Min. :14.89
## 1st Qu.:102.7 1st Qu.: 9.06 1st Qu.:8.12 1st Qu.:26.76
## Median :113.1 Median :10.02 Median :8.24 Median :29.14
## Mean :113.3 Mean :10.16 Mean :8.26 Mean :28.50
## 3rd Qu.:124.3 3rd Qu.:11.20 3rd Qu.:8.42 3rd Qu.:30.83
## Max. :177.3 Max. :16.85 Max. :9.20 Max. :33.88
##
## J.date SIN_TIME COS_TIME Round
## Min. : 79.0 Min. :-0.9882 Min. :-0.9953 Min. : 1.000
## 1st Qu.:134.0 1st Qu.:-0.1606 1st Qu.:-0.9765 1st Qu.: 5.000
## Median :153.0 Median : 0.4878 Median :-0.8278 Median : 7.000
## Mean :162.7 Mean : 0.3044 Mean :-0.6873 Mean : 6.271
## 3rd Qu.:192.0 3rd Qu.: 0.7423 3rd Qu.:-0.5037 3rd Qu.: 8.000
## Max. :286.0 Max. : 1.0000 Max. : 0.2102 Max. :10.000
##
## Habitat Site Set
## Eelgrass :1435 E5 :275 Min. :1.000
## Sand flat: 0 E3 :265 1st Qu.:1.000
## E2 :255 Median :2.000
## E4 :235 Mean :2.019
## E1 :205 3rd Qu.:3.000
## E7 :105 Max. :4.000
## (Other): 95 NA's :1
## Species Life_stage abundance
## Three-spined stickleback:245 0 : 17 Min. : 0.00
## Shiner surfperch :242 adult :665 1st Qu.: 1.00
## Starry flounder :162 juvenile:539 Median : 2.00
## Bay pipefish :131 NA's :214 Mean : 26.05
## Tubesnout : 83 3rd Qu.: 10.00
## Surf smelt : 74 Max. :1698.00
## (Other) :498
## class round month
## 0 : 17 Min. : 1.000 April :215
## migratory: 214 1st Qu.: 5.000 July :231
## resident :1185 Median : 6.000 June :347
## NA's : 19 Mean : 6.147 March : 95
## 3rd Qu.: 8.000 May :418
## Max. :10.000 October : 52
## September: 77sapply(eelgrass, sd)## Year Date Temp.surf DO.surf DOmg.surf pH.surf
## 0.4980858 13.9351400 2.3601026 15.3515577 1.3885806 0.2706839
## Sal.surf Temp.bot DO.bot DOmg.bot pH.bot Sal.bot
## 6.0302167 2.2570560 18.3091380 1.7728437 0.2680020 3.4087278
## J.date SIN_TIME COS_TIME Round Habitat Site
## 45.5103676 0.5526942 0.3602235 2.2162626 0.0000000 1.7656842
## Set Species Life_stage abundance class round
## NA 11.6310130 NA 95.3126364 NA 2.0735395
## month
## 1.7078374# Sand flat summary
sandflat <- subset(purse, Habitat %in% "Sand flat")
summary(sandflat)## Year Date Temp.surf DO.surf
## Min. :0.0000 6/3/2017 : 90 Min. : 7.55 Min. : 43.1
## 1st Qu.:0.0000 7/14/2017: 57 1st Qu.:12.62 1st Qu.: 98.4
## Median :0.0000 12-Jun-16: 51 Median :14.21 Median :109.9
## Mean :0.4369 6/19/2017: 50 Mean :14.38 Mean :106.8
## 3rd Qu.:1.0000 5/18/2017: 48 3rd Qu.:15.75 3rd Qu.:116.0
## Max. :1.0000 5/6/2017 : 42 Max. :20.41 Max. :143.5
## (Other) :454
## DOmg.surf pH.surf Sal.surf Temp.bot
## Min. : 4.04 Min. :7.810 Min. : 0.55 Min. : 6.19
## 1st Qu.: 9.41 1st Qu.:8.015 1st Qu.: 5.27 1st Qu.:12.35
## Median :10.36 Median :8.140 Median :10.97 Median :13.65
## Mean :10.28 Mean :8.187 Mean :11.70 Mean :13.83
## 3rd Qu.:11.44 3rd Qu.:8.310 3rd Qu.:16.42 3rd Qu.:15.50
## Max. :14.56 Max. :8.700 Max. :30.84 Max. :20.22
##
## DO.bot DOmg.bot pH.bot Sal.bot
## Min. : 45.3 Min. : 4.130 Min. :7.770 Min. : 1.31
## 1st Qu.: 97.9 1st Qu.: 9.258 1st Qu.:8.020 1st Qu.:12.17
## Median :108.2 Median : 9.740 Median :8.120 Median :16.14
## Mean :107.1 Mean : 9.968 Mean :8.174 Mean :17.55
## 3rd Qu.:114.7 3rd Qu.:11.140 3rd Qu.:8.280 3rd Qu.:24.35
## Max. :144.2 Max. :13.360 Max. :8.670 Max. :34.05
##
## J.date SIN_TIME COS_TIME Round
## Min. : 79.0 Min. :-0.99771 Min. :-0.9999 Min. : 1.000
## 1st Qu.:126.0 1st Qu.: 0.01075 1st Qu.:-0.9696 1st Qu.: 4.000
## Median :154.0 Median : 0.47276 Median :-0.8644 Median : 6.000
## Mean :158.7 Mean : 0.35872 Mean :-0.6889 Mean : 5.996
## 3rd Qu.:182.0 3rd Qu.: 0.82719 3rd Qu.:-0.5185 3rd Qu.: 8.000
## Max. :285.0 Max. : 0.99999 Max. : 0.2102 Max. :10.000
##
## Habitat Site Set
## Eelgrass : 0 SF6 :178 Min. :1.000
## Sand flat:792 SF1 :137 1st Qu.:1.000
## SF3 :135 Median :2.000
## SF4 :121 Mean :1.986
## SF2 :117 3rd Qu.:3.000
## SF5 :104 Max. :3.000
## (Other): 0
## Species Life_stage abundance
## Starry flounder :131 0 : 65 Min. : 0.000
## Shiner surfperch : 87 adult :219 1st Qu.: 1.000
## Three-spined stickleback: 75 juvenile:423 Median : 1.000
## Chinook : 74 NA's : 85 Mean : 8.374
## 0 : 65 3rd Qu.: 3.250
## Arrow goby : 64 Max. :626.000
## (Other) :296
## class round month
## 0 : 65 Min. : 1.000 April :142
## migratory:156 1st Qu.: 4.000 July :110
## resident :532 Median : 6.000 June :191
## NA's : 39 Mean : 5.951 March : 58
## 3rd Qu.: 7.500 May :238
## Max. :10.000 October : 28
## September: 25sapply(sandflat, sd)## Year Date Temp.surf DO.surf DOmg.surf pH.surf
## 0.4963119 14.3669995 2.7012171 14.6284763 1.5976603 0.2170129
## Sal.surf Temp.bot DO.bot DOmg.bot pH.bot Sal.bot
## 7.6917446 2.6065195 15.1503214 1.4979394 0.2192653 7.9716065
## J.date SIN_TIME COS_TIME Round Habitat Site
## 43.2564944 0.5225303 0.3523348 2.2218854 0.0000000 1.7858508
## Set Species Life_stage abundance class round
## 0.8083377 13.6352066 NA 37.4836395 NA 2.0455862
## month
## 1.6705429### Grab just chinook catch - automatically removes any 0s or unidentified
### species. 1
p.1 <- purse[which(purse$Species == "Chinook"), ]
### standardize variables install.packages('robustHD')
### install.packages('robustbase', repos='http://R-Forge.R-project.org')
library(robustHD)
p.1$s.temp <- standardize(p.1$Temp.surf, centerFun = mean, scaleFun = sd)
summary(p.1$s.temp) #can see that mean is now 0, and SD is on the original scale of x (temp) -- i.e. predictor centering ## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -2.9418 -0.5635 -0.1433 0.0000 0.4847 2.8797p.1$s.sal <- standardize(p.1$Sal.surf, centerFun = mean, scaleFun = sd)
p.1$s.do <- standardize(p.1$DOmg.surf, centerFun = mean, scaleFun = sd)
p.1$s.pH <- standardize(p.1$pH.surf, centerFun = mean, scaleFun = sd)
p.1$s.J.date <- standardize(p.1$J.date, centerFun = mean, scaleFun = sd)
### Create variable j2 which is Julian day squared in order to represent J
### date as quadratic relationship instead of linear
p.1$j2 <- p.1$s.J.date^2
# summarize by site-day
library(plyr)
p.chin <- ddply(p.1, .(Year, J.date, s.J.date, j2, Habitat, Site, s.temp, s.sal,
s.do, s.pH), summarize, abundance = sum(abundance))
# plot abundance~J.date to see data distribution
plot(abundance ~ J.date, data = p.chin)
Habitat variables
Add eelgrass and sand flat habitat variables for each site
p.hab <- read.csv("/Users/Lia/Documents/Git/Fraser-salmon/all.data/site.char.eelgrass.csv")
summary(p.hab)## Collection.date Site.ID Site shoot_den
## :7 :6 E1 :1 Min. : 0.0
## 18-Jul-16 :2 017 FR :1 E2 :1 1st Qu.: 0.0
## 20-Jul-16 :1 017BFR :1 E3 :1 Median : 304.0
## archipelago:3 018 FR :1 E4 :1 Mean : 349.6
## 030 FR :1 E5 :1 3rd Qu.: 508.4
## WP126(19 FR):1 E6 :1 Max. :1118.9
## (Other) :2 (Other):7
## percent_cov leaf_area_index epiphyte_wt meanturb
## Min. : 0.000 Min. :0.000 Min. : 0.000 Min. : 2.017
## 1st Qu.: 0.000 1st Qu.:0.000 1st Qu.: 0.000 1st Qu.: 3.465
## Median : 7.667 Median :2.000 Median : 0.000 Median : 8.971
## Mean : 39.863 Mean :2.017 Mean : 1.788 Mean : 8.113
## 3rd Qu.: 75.000 3rd Qu.:3.451 3rd Qu.: 0.484 3rd Qu.:10.941
## Max. :100.000 Max. :5.831 Max. :10.112 Max. :17.810
## NA's :3sapply(p.hab, sd)## Warning in var(if (is.vector(x) || is.factor(x)) x else as.double(x), na.rm = na.rm): Calling var(x) on a factor x is deprecated and will become an error.
## Use something like 'all(duplicated(x)[-1L])' to test for a constant vector.
## Warning in var(if (is.vector(x) || is.factor(x)) x else as.double(x), na.rm = na.rm): Calling var(x) on a factor x is deprecated and will become an error.
## Use something like 'all(duplicated(x)[-1L])' to test for a constant vector.
## Warning in var(if (is.vector(x) || is.factor(x)) x else as.double(x), na.rm = na.rm): Calling var(x) on a factor x is deprecated and will become an error.
## Use something like 'all(duplicated(x)[-1L])' to test for a constant vector.## Collection.date Site.ID Site shoot_den
## 1.290994 2.577019 3.894440 415.152611
## percent_cov leaf_area_index epiphyte_wt meanturb
## 44.244763 2.217159 NA 5.011330# eelgrass summary
eelgrass <- subset(p.hab, Site %in% c("E1", "E2", "E3", "E4", "E5", "E6", "E7"))
summary(eelgrass)## Collection.date Site.ID Site shoot_den
## :1 017 FR :1 E1 :1 Min. : 304.0
## 18-Jul-16 :2 017BFR :1 E2 :1 1st Qu.: 400.0
## 20-Jul-16 :1 018 FR :1 E3 :1 Median : 508.4
## archipelago:3 030 FR :1 E4 :1 Mean : 649.2
## WP126(19 FR):1 E5 :1 3rd Qu.: 906.6
## WP227(20 FR):1 E6 :1 Max. :1118.9
## (Other) :1 (Other):1
## percent_cov leaf_area_index epiphyte_wt meanturb
## Min. : 7.667 Min. :2.000 Min. : 0.484 Min. :2.017
## 1st Qu.: 75.000 1st Qu.:2.815 1st Qu.: 0.484 1st Qu.:2.295
## Median : 75.000 Median :3.451 Median : 3.640 Median :3.465
## Mean : 74.032 Mean :3.746 Mean : 4.469 Mean :4.430
## 3rd Qu.: 92.778 3rd Qu.:4.656 3rd Qu.: 7.625 3rd Qu.:5.820
## Max. :100.000 Max. :5.831 Max. :10.112 Max. :9.295
## NA's :3sapply(eelgrass, sd)## Warning in var(if (is.vector(x) || is.factor(x)) x else as.double(x), na.rm = na.rm): Calling var(x) on a factor x is deprecated and will become an error.
## Use something like 'all(duplicated(x)[-1L])' to test for a constant vector.
## Warning in var(if (is.vector(x) || is.factor(x)) x else as.double(x), na.rm = na.rm): Calling var(x) on a factor x is deprecated and will become an error.
## Use something like 'all(duplicated(x)[-1L])' to test for a constant vector.
## Warning in var(if (is.vector(x) || is.factor(x)) x else as.double(x), na.rm = na.rm): Calling var(x) on a factor x is deprecated and will become an error.
## Use something like 'all(duplicated(x)[-1L])' to test for a constant vector.## Collection.date Site.ID Site shoot_den
## 1.214986 2.160247 2.160247 343.146423
## percent_cov leaf_area_index epiphyte_wt meanturb
## 31.049192 1.508139 NA 2.771784# Sand flat summary
sandflat <- subset(p.hab, Site %in% c("SF1", "SF2", "SF3", "SF4", "SF5", "SF6"))
summary(sandflat)## Collection.date Site.ID Site shoot_den percent_cov
## :6 :6 SF1 :1 Min. :0 Min. :0
## 18-Jul-16 :0 017 FR :0 SF2 :1 1st Qu.:0 1st Qu.:0
## 20-Jul-16 :0 017BFR :0 SF3 :1 Median :0 Median :0
## archipelago:0 018 FR :0 SF4 :1 Mean :0 Mean :0
## 030 FR :0 SF5 :1 3rd Qu.:0 3rd Qu.:0
## WP126(19 FR):0 SF6 :1 Max. :0 Max. :0
## (Other) :0 (Other):0
## leaf_area_index epiphyte_wt meanturb
## Min. :0 Min. :0 Min. : 8.971
## 1st Qu.:0 1st Qu.:0 1st Qu.:10.427
## Median :0 Median :0 Median :11.901
## Mean :0 Mean :0 Mean :12.411
## 3rd Qu.:0 3rd Qu.:0 3rd Qu.:13.435
## Max. :0 Max. :0 Max. :17.810
## sapply(sandflat, sd)## Warning in var(if (is.vector(x) || is.factor(x)) x else as.double(x), na.rm = na.rm): Calling var(x) on a factor x is deprecated and will become an error.
## Use something like 'all(duplicated(x)[-1L])' to test for a constant vector.
## Warning in var(if (is.vector(x) || is.factor(x)) x else as.double(x), na.rm = na.rm): Calling var(x) on a factor x is deprecated and will become an error.
## Use something like 'all(duplicated(x)[-1L])' to test for a constant vector.
## Warning in var(if (is.vector(x) || is.factor(x)) x else as.double(x), na.rm = na.rm): Calling var(x) on a factor x is deprecated and will become an error.
## Use something like 'all(duplicated(x)[-1L])' to test for a constant vector.## Collection.date Site.ID Site shoot_den
## 0.000000 0.000000 1.870829 0.000000
## percent_cov leaf_area_index epiphyte_wt meanturb
## 0.000000 0.000000 0.000000 3.145279## Now add to p.chin
p.hab2 <- p.hab[, c(3:6, 8)]
p.chin.hab <- p.chin
p.chin.hab$leaf_area_index <- p.hab2[match(p.chin.hab$Site, p.hab2$Site), 4]
p.chin.hab$meanturb <- p.hab2[match(p.chin.hab$Site, p.hab2$Site), 5]
## standardize
p.chin.hab$leaf_area_index <- standardize(p.chin.hab$leaf_area_index, centerFun = mean,
scaleFun = sd)
p.chin.hab$leaf_area_index <- as.numeric(p.chin.hab$leaf_area_index)
p.chin.hab$meanturb <- standardize(p.chin.hab$meanturb, centerFun = mean, scaleFun = sd)
p.chin.hab$meanturb <- as.numeric(p.chin.hab$meanturb)VIF for collinearity of habitat variables
Assess variance inflation factors
library(car)
library(GGally)
vif(lm(abundance ~ Habitat + leaf_area_index + s.J.date + j2 + meanturb + Year +
s.temp + s.sal + s.do + s.pH, data = p.chin.hab))
## Habitat leaf_area_index s.J.date j2
## 7.160505 6.972103 3.050891 2.062968
## meanturb Year s.temp s.sal
## 4.898655 1.358753 3.174750 5.458973
## s.do s.pH
## 2.720738 1.968218
# alias function identifies covariates that are multiples of each other - in
# this case we are OK.
alias(lm(abundance ~ Habitat + leaf_area_index + s.J.date + j2 + meanturb +
Year + s.temp + s.sal + s.do + s.pH, data = p.chin.hab))
## Model :
## abundance ~ Habitat + leaf_area_index + s.J.date + j2 + meanturb +
## Year + s.temp + s.sal + s.do + s.pH
## Pearson Corr with all vars
year <- p.chin.hab$Year
jday <- p.chin.hab$s.J.date
j2 <- p.chin.hab$j2
temp <- p.chin.hab$s.temp
do <- p.chin.hab$s.do
ph <- p.chin.hab$s.pH
sal <- p.chin.hab$s.sal
lai <- p.chin.hab$leaf_area_index
turb <- p.chin.hab$meanturb
hab <- p.chin.hab$Habitat
habcovar <- cbind.data.frame(hab, year, lai, turb, jday, j2, temp, do, ph, sal)
ggpairs(data = na.omit(habcovar), title = "Pearson Correlation plot habitat variables")
ggsave("/Users/Lia/Documents/Git/Fraser-salmon/Ch1.Seasonal_diversity/expl.plots/PearsonCorrAllVars_chinook.pdf")Model selection
Full model
We started with a full global model with all the abiotic factors and all the habitat covariates included as well as site as a random effect. That model did not converge, so we performed manual backward selection beginning with the variables with highest VIF, and retaining the top variables from the corresponding marsh model, until the model converged. Throughout the process we viewed the VIF and model assumption plots, and used a hypothesis-based approach to select variables (as well as trying multiple iterations).
library(MASS)
library(lme4)## Warning: package 'lme4' was built under R version 3.4.4## Chinook: full model with habitat variables
p.chin1 <- glmer.nb(abundance ~ Habitat + s.J.date + s.temp + Year + s.do +
s.pH + (1 | Site), data = p.chin.hab, na.action = "na.fail")
summary(p.chin1) ## AIC 365.7 ## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: Negative Binomial(4.4747) ( log )
## Formula: abundance ~ Habitat + s.J.date + s.temp + Year + s.do + s.pH +
## (1 | Site)
## Data: p.chin.hab
##
## AIC BIC logLik deviance df.resid
## 365.7 387.8 -173.9 347.7 77
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.4021 -0.5944 -0.2815 0.4967 3.1490
##
## Random effects:
## Groups Name Variance Std.Dev.
## Site (Intercept) 0.1119 0.3345
## Number of obs: 86, groups: Site, 13
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.98589 0.23431 4.208 2.58e-05 ***
## HabitatSand flat -0.51129 0.29608 -1.727 0.08419 .
## s.J.date -0.12703 0.13666 -0.930 0.35261
## s.temp 0.07661 0.12192 0.628 0.52977
## Year 0.52378 0.18808 2.785 0.00535 **
## s.do 0.26861 0.11149 2.409 0.01598 *
## s.pH -0.03519 0.12091 -0.291 0.77100
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) HbttSf s.J.dt s.temp Year s.do
## HabttSndflt -0.680
## s.J.date -0.034 0.068
## s.temp 0.185 -0.286 -0.604
## Year -0.533 0.084 -0.057 0.121
## s.do 0.027 -0.112 0.459 -0.084 0.103
## s.pH -0.299 0.467 -0.009 -0.331 -0.051 -0.182## Introduce function to calculate Variance Inflation Factors for GLMMs
vif.lme <- function(fit) {
## adapted from rms::vif
v <- vcov(fit)
nam <- names(fixef(fit))
## exclude intercepts
ns <- sum(1 * (nam == "Intercept" | nam == "(Intercept)"))
if (ns > 0) {
v <- v[-(1:ns), -(1:ns), drop = FALSE]
nam <- nam[-(1:ns)]
}
d <- diag(v)^0.5
v <- diag(solve(v/(d %o% d)))
names(v) <- nam
v
}
# test VIF of full model -- mostly OK
vif.lme(p.chin1)## HabitatSand flat s.J.date s.temp Year
## 1.345749 2.229278 2.053339 1.053817
## s.do s.pH
## 1.409247 1.449275# Test fit and assumptions of full (global) model Plot residuals vs. fitted
# values
plot(fitted(p.chin1), resid(p.chin1), main = "Global Chinook GLMM", xlab = "Fitted",
ylab = "Pearson residuals")
## q plot
qqnorm(x = p.chin.hab$abundance, y = resid(p.chin1), main = "Q-Q Global Chinook GLMM")
qqline(resid(p.chin1), col = 2)
library(sjPlot)## Warning: package 'sjPlot' was built under R version 3.4.3## Warning in checkMatrixPackageVersion(): Package version inconsistency detected.
## TMB was built with Matrix version 1.2.12
## Current Matrix version is 1.2.11
## Please re-install 'TMB' from source using install.packages('TMB', type = 'source') or ask CRAN for a binary version of 'TMB' matching CRAN's 'Matrix' package# QQ plot of random effects quantiles against standard normal quantiles
sjp.glmer(p.chin1, type = "re.qq")
# Check model assumptions - note that we are not actually assuming linear
# relationships with Jday, but can't tell this general function that.
sjp.glmer(p.chin1, type = "ma")
# See marginal and conditional pseudo R2 for full model
library(piecewiseSEM)
rsquared(p.chin1, aicc = TRUE)## Class Family Link n Marginal Conditional
## 1 glmerMod Negative Binomial log 86 0.3587852 0.5289181Model selection - dredge
Then continued to dredge function to determine optimal model using AICc (chin has somewhat small sample size: 86 observations at 13 sites).
library(MuMIn)
library(knitr)
# Generate model set
p.model.set.chin <- dredge(p.chin1)
# p.model.set.chin.AIC<- dredge(p.chin1, rank = 'AIC') #compare to AIC to
# see if much different
# Create model selection table
p.model_table.chin <- model.sel(p.model.set.chin)
options(scipen = 7)
names(p.model_table.chin) <- c("(Int)", "Habitat", "DO", "Jday", "pH", "Temp",
"Year", "df", "LL", "AICc", "delta", "weight")
kable(head(p.model_table.chin, n = 100), digits = 3)| (Int) | Habitat | DO | Jday | pH | Temp | Year | df | LL | AICc | delta | weight | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 36 | 0.956 | + | 0.313 | NA | NA | NA | 0.514 | 6 | -174.344 | 361.751 | 0.000 | 0.20456373180 |
| 35 | 0.705 | NA | 0.308 | NA | NA | NA | 0.544 | 5 | -175.782 | 362.314 | 0.563 | 0.15440510347 |
| 40 | 0.953 | + | 0.273 | -0.078 | NA | NA | 0.510 | 7 | -174.060 | 363.556 | 1.805 | 0.08297937494 |
| 39 | 0.705 | NA | 0.262 | -0.087 | NA | NA | 0.539 | 6 | -175.433 | 363.929 | 2.177 | 0.06886675942 |
| 44 | 0.972 | + | 0.316 | NA | -0.032 | NA | 0.514 | 7 | -174.302 | 364.039 | 2.288 | 0.06516961517 |
| 52 | 0.956 | + | 0.312 | NA | NA | -0.003 | 0.513 | 7 | -174.343 | 364.123 | 2.371 | 0.06250457753 |
| 51 | 0.705 | NA | 0.302 | NA | NA | -0.026 | 0.539 | 6 | -175.739 | 364.542 | 2.790 | 0.05068528339 |
| 43 | 0.705 | NA | 0.307 | NA | 0.019 | NA | 0.544 | 6 | -175.768 | 364.599 | 2.847 | 0.04926704491 |
| 56 | 0.965 | + | 0.263 | -0.127 | NA | 0.065 | 0.521 | 8 | -173.899 | 365.669 | 3.917 | 0.02885127572 |
| 48 | 0.958 | + | 0.275 | -0.076 | -0.010 | NA | 0.510 | 8 | -174.056 | 365.982 | 4.231 | 0.02466810618 |
| 47 | 0.705 | NA | 0.253 | -0.098 | 0.045 | NA | 0.539 | 7 | -175.354 | 366.145 | 4.393 | 0.02274471869 |
| 55 | 0.705 | NA | 0.257 | -0.113 | NA | 0.034 | 0.545 | 7 | -175.388 | 366.212 | 4.461 | 0.02199037633 |
| 60 | 0.977 | + | 0.318 | NA | -0.036 | 0.009 | 0.516 | 8 | -174.297 | 366.465 | 4.713 | 0.01938115146 |
| 59 | 0.705 | NA | 0.297 | NA | 0.037 | -0.037 | 0.536 | 7 | -175.691 | 366.818 | 5.067 | 0.01624125323 |
| 4 | 1.310 | + | 0.288 | NA | NA | NA | NA | 5 | -178.186 | 367.122 | 5.371 | 0.01395158938 |
| 37 | 0.726 | NA | NA | -0.213 | NA | NA | 0.489 | 5 | -178.409 | 367.569 | 5.817 | 0.01115858447 |
| 38 | 0.955 | + | NA | -0.208 | NA | NA | 0.461 | 6 | -177.270 | 367.604 | 5.852 | 0.01096561749 |
| 64 | 0.986 | + | 0.269 | -0.127 | -0.035 | 0.077 | 0.524 | 9 | -173.857 | 368.082 | 6.331 | 0.00863202856 |
| 3 | 1.022 | NA | 0.281 | NA | NA | NA | NA | 4 | -179.981 | 368.457 | 6.705 | 0.00715789797 |
| 63 | 0.705 | NA | 0.251 | -0.113 | 0.038 | 0.022 | 0.543 | 8 | -175.337 | 368.545 | 6.793 | 0.00685048720 |
| 8 | 1.305 | + | 0.242 | -0.088 | NA | NA | NA | 6 | -177.819 | 368.701 | 6.950 | 0.00633389601 |
| 54 | 0.976 | + | NA | -0.281 | NA | 0.109 | 0.482 | 7 | -176.789 | 369.013 | 7.262 | 0.00541882828 |
| 45 | 0.725 | NA | NA | -0.226 | 0.098 | NA | 0.490 | 6 | -178.010 | 369.084 | 7.333 | 0.00523102988 |
| 20 | 1.301 | + | 0.283 | NA | NA | -0.027 | NA | 6 | -178.140 | 369.344 | 7.592 | 0.00459433372 |
| 53 | 0.725 | NA | NA | -0.266 | NA | 0.080 | 0.505 | 6 | -178.142 | 369.346 | 7.595 | 0.00458792518 |
| 12 | 1.326 | + | 0.292 | NA | -0.031 | NA | NA | 6 | -178.144 | 369.351 | 7.599 | 0.00457835700 |
| 46 | 0.927 | + | NA | -0.216 | 0.055 | NA | 0.461 | 7 | -177.147 | 369.730 | 7.978 | 0.00378748302 |
| 7 | 1.020 | NA | 0.230 | -0.095 | NA | NA | NA | 5 | -179.557 | 369.864 | 8.112 | 0.00354241144 |
| 19 | 1.017 | NA | 0.271 | NA | NA | -0.047 | NA | 5 | -179.840 | 370.429 | 8.678 | 0.00266946629 |
| 11 | 1.022 | NA | 0.279 | NA | 0.013 | NA | NA | 5 | -179.974 | 370.698 | 8.947 | 0.00233361241 |
| 24 | 1.314 | + | 0.237 | -0.111 | NA | 0.032 | NA | 7 | -177.779 | 370.994 | 9.243 | 0.00201258153 |
| 16 | 1.308 | + | 0.244 | -0.086 | -0.007 | NA | NA | 7 | -177.817 | 371.070 | 9.318 | 0.00193809015 |
| 34 | 0.965 | + | NA | NA | NA | NA | 0.457 | 5 | -180.173 | 371.096 | 9.344 | 0.00191330723 |
| 61 | 0.724 | NA | NA | -0.258 | 0.080 | 0.052 | 0.501 | 7 | -177.910 | 371.256 | 9.505 | 0.00176582713 |
| 33 | 0.731 | NA | NA | NA | NA | NA | 0.482 | 4 | -181.446 | 371.385 | 9.634 | 0.00165528238 |
| 62 | 0.966 | + | NA | -0.280 | 0.017 | 0.102 | 0.481 | 8 | -176.779 | 371.427 | 9.676 | 0.00162075377 |
| 6 | 1.275 | + | NA | -0.205 | NA | NA | NA | 5 | -180.404 | 371.557 | 9.806 | 0.00151883333 |
| 28 | 1.314 | + | 0.287 | NA | -0.021 | -0.019 | NA | 7 | -178.124 | 371.684 | 9.933 | 0.00142544201 |
| 15 | 1.020 | NA | 0.220 | -0.105 | 0.042 | NA | NA | 6 | -179.488 | 372.039 | 10.288 | 0.00119358174 |
| 23 | 1.021 | NA | 0.229 | -0.098 | NA | 0.004 | NA | 6 | -179.556 | 372.176 | 10.424 | 0.00111490008 |
| 5 | 1.013 | NA | NA | -0.208 | NA | NA | NA | 4 | -181.879 | 372.251 | 10.499 | 0.00107376227 |
| 27 | 1.015 | NA | 0.264 | NA | 0.046 | -0.063 | NA | 6 | -179.765 | 372.594 | 10.842 | 0.00090454585 |
| 49 | 0.728 | NA | NA | NA | NA | -0.078 | 0.471 | 5 | -181.003 | 372.756 | 11.004 | 0.00083425616 |
| 50 | 0.949 | + | NA | NA | NA | -0.061 | 0.448 | 6 | -179.912 | 372.887 | 11.136 | 0.00078114591 |
| 22 | 1.300 | + | NA | -0.253 | NA | 0.076 | NA | 6 | -180.158 | 373.379 | 11.628 | 0.00061073715 |
| 41 | 0.732 | NA | NA | NA | 0.054 | NA | 0.482 | 5 | -181.318 | 373.386 | 11.635 | 0.00060864901 |
| 42 | 0.960 | + | NA | NA | 0.011 | NA | 0.457 | 6 | -180.167 | 373.398 | 11.647 | 0.00060498902 |
| 32 | 1.327 | + | 0.240 | -0.111 | -0.021 | 0.039 | NA | 8 | -177.764 | 373.399 | 11.648 | 0.00060478186 |
| 14 | 1.250 | + | NA | -0.212 | 0.052 | NA | NA | 6 | -180.287 | 373.637 | 11.886 | 0.00053689177 |
| 13 | 1.014 | NA | NA | -0.220 | 0.093 | NA | NA | 5 | -181.519 | 373.789 | 12.037 | 0.00049766863 |
| 57 | 0.730 | NA | NA | NA | 0.105 | -0.109 | 0.464 | 6 | -180.569 | 374.200 | 12.449 | 0.00040510387 |
| 21 | 1.019 | NA | NA | -0.238 | NA | 0.049 | NA | 5 | -181.774 | 374.298 | 12.547 | 0.00038577707 |
| 31 | 1.019 | NA | 0.222 | -0.098 | 0.046 | -0.012 | NA | 7 | -179.483 | 374.401 | 12.650 | 0.00036636411 |
| 2 | 1.278 | + | NA | NA | NA | NA | NA | 4 | -183.286 | 375.066 | 13.314 | 0.00026280648 |
| 58 | 0.921 | + | NA | NA | 0.050 | -0.077 | 0.445 | 7 | -179.818 | 375.072 | 13.320 | 0.00026202510 |
| 30 | 1.284 | + | NA | -0.251 | 0.027 | 0.066 | NA | 7 | -180.131 | 375.699 | 13.947 | 0.00019152803 |
| 1 | 1.013 | NA | NA | NA | NA | NA | NA | 3 | -184.837 | 375.966 | 14.214 | 0.00016757094 |
| 29 | 1.016 | NA | NA | -0.229 | 0.086 | 0.016 | NA | 6 | -181.510 | 376.082 | 14.331 | 0.00015809118 |
| 18 | 1.253 | + | NA | NA | NA | -0.075 | NA | 5 | -182.893 | 376.536 | 14.784 | 0.00012602500 |
| 17 | 1.003 | NA | NA | NA | NA | -0.092 | NA | 4 | -184.246 | 376.986 | 15.235 | 0.00010059887 |
| 10 | 1.273 | + | NA | NA | 0.011 | NA | NA | 5 | -183.281 | 377.312 | 15.560 | 0.00008549940 |
| 9 | 1.013 | NA | NA | NA | 0.049 | NA | NA | 4 | -184.730 | 377.953 | 16.201 | 0.00006204732 |
| 25 | 1.000 | NA | NA | NA | 0.113 | -0.127 | NA | 5 | -183.753 | 378.256 | 16.504 | 0.00005332741 |
| 26 | 1.218 | + | NA | NA | 0.060 | -0.094 | NA | 6 | -182.760 | 378.583 | 16.831 | 0.00004528767 |
# determine at which delta AICc score you reach cumulative AICc weight of
# 0.95 (Harrison et al. 2018)
find_cumsum = function(df, delta) {
for (i in nrow(df)) {
df = df[df$delta >= delta, ]
return(min(df$delta[cumsum(df$weight) >= 0.95]))
}
}
find_cumsum(p.model_table.chin, 0) ##THIS SUGGESTS AVERAGING UP TO delta AICc 7.59. We will use a cutoff of delta AICc < 4.
## [1] 7.592112
# Model averaging Version 1: use all models with delta AIC score less than 4
p.model.set.chin.4 <- get.models(p.model.set.chin, subset = delta < 4)
p.avg_model.chin.4 <- model.avg(p.model.set.chin.4)
summary(p.avg_model.chin.4)
##
## Call:
## model.avg(object = p.model.set.chin.4)
##
## Component model call:
## glmer(formula = abundance ~ <9 unique rhs>, data = p.chin.hab,
## family = negative.binomial(theta = 4.47467009362264), na.action =
## na.fail)
##
## Component models:
## df logLik AICc delta weight
## 126 6 -174.34 361.75 0.00 0.27
## 26 5 -175.78 362.31 0.56 0.20
## 1236 7 -174.06 363.56 1.80 0.11
## 236 6 -175.43 363.93 2.18 0.09
## 1246 7 -174.30 364.04 2.29 0.08
## 1256 7 -174.34 364.12 2.37 0.08
## 256 6 -175.74 364.54 2.79 0.07
## 246 6 -175.77 364.60 2.85 0.06
## 12356 8 -173.90 365.67 3.92 0.04
##
## Term codes:
## Habitat s.do s.J.date s.pH s.temp Year
## 1 2 3 4 5 6
##
## Model-averaged coefficients:
## (full average)
## Estimate Std. Error Adjusted SE z value Pr(>|z|)
## (Intercept) 0.8516625 0.2406532 0.2433611 3.500 0.000466 ***
## HabitatSand flat -0.2651547 0.2999038 0.3019090 0.878 0.379802
## s.do 0.3002695 0.0993223 0.1007970 2.979 0.002892 **
## Year 0.5255259 0.1874715 0.1903397 2.761 0.005763 **
## s.J.date -0.0210789 0.0659719 0.0666429 0.316 0.751778
## s.pH -0.0015055 0.0434506 0.0440854 0.034 0.972758
## s.temp 0.0004921 0.0428127 0.0434058 0.011 0.990954
##
## (conditional average)
## Estimate Std. Error Adjusted SE z value Pr(>|z|)
## (Intercept) 0.851662 0.240653 0.243361 3.500 0.000466 ***
## HabitatSand flat -0.458153 0.258816 0.262814 1.743 0.081288 .
## s.do 0.300269 0.099322 0.100797 2.979 0.002892 **
## Year 0.525526 0.187472 0.190340 2.761 0.005763 **
## s.J.date -0.089507 0.111159 0.112846 0.793 0.427672
## s.pH -0.010094 0.112125 0.113774 0.089 0.929303
## s.temp 0.002658 0.099476 0.100855 0.026 0.978971
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Relative variable importance:
## s.do Year Habitat s.J.date s.temp s.pH
## Importance: 1.00 1.00 0.58 0.24 0.19 0.15
## N containing models: 9 9 5 3 3 2
p.chin.ci <- data.frame(confint(p.avg_model.chin.4, full = TRUE))
# Get pseudo R squared values for models up to delta < 4 model.list.chin<-
# list(#manually list top x models from dredge)
p.chin.4.Rsq <- rsquared(p.model.set.chin.4, aicc = TRUE)
## write tables to .csv for easy comparison and plugging into appendix table
p.avg_model_4df.chin <- data.frame(p.avg_model.chin.4$msTable)
p.avg_model_components4.chin <- cbind(p.chin.4.Rsq, p.avg_model_4df.chin)
p.r = data.frame(Coeff = rownames(p.avg_model_4df.chin, rep(NA, length(p.avg_model_components4.chin))))
p.avg_model_components4.chin <- cbind(p.avg_model_components4.chin, p.r)
p.avg_model_components4.chin <- p.avg_model_components4.chin[, -c(7, 8)]
# write.csv(p.avg_model_components4.chin,
# '/Users/Lia/Documents/Git/Fraser-salmon/all.data/avg_model_components4_pursechin.csv')Parameter Plot
The results of model averaging including all top ranked models up to delta AICc 4
## Warning: package 'cowplot' was built under R version 3.4.3##
## Attaching package: 'cowplot'## The following objects are masked from 'package:sjPlot':
##
## plot_grid, save_plot## The following object is masked from 'package:ggplot2':
##
## ggsave
Check fit of top model by avg – including all parameters with weight >0.5
# Test fit and assumptions of final (averaged) model, following the nesting
# rule
p.avg.chin <- glmer.nb(abundance ~ s.do + Year + Habitat + (1 | Site), data = p.chin.hab,
na.action = "na.fail")
summary(p.avg.chin) # AIC 360.7## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: Negative Binomial(4.3106) ( log )
## Formula: abundance ~ s.do + Year + Habitat + (1 | Site)
## Data: p.chin.hab
##
## AIC BIC logLik deviance df.resid
## 360.7 375.4 -174.3 348.7 80
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.3717 -0.6216 -0.2597 0.3799 3.1145
##
## Random effects:
## Groups Name Variance Std.Dev.
## Site (Intercept) 0.1021 0.3195
## Number of obs: 86, groups: Site, 13
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.95601 0.21976 4.350 0.0000136 ***
## s.do 0.31274 0.09416 3.322 0.000895 ***
## Year 0.51413 0.18717 2.747 0.006017 **
## HabitatSand flat -0.45215 0.25286 -1.788 0.073757 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) s.do Year
## s.do -0.032
## Year -0.602 0.112
## HabttSndflt -0.634 -0.047 0.148vif.lme(p.avg.chin)## s.do Year HabitatSand flat
## 1.016995 1.037407 1.026673### Plot residuals vs. fitted values
plot(fitted(p.avg.chin), resid(p.avg.chin), main = "Averaged Chinook GLMM",
xlab = "Fitted", ylab = "Pearson residuals")
## q plot
qqnorm(x = p.chin.hab$abundance, y = resid(p.avg.chin), main = "Q-Q Averaged Chinook GLMM")
qqline(resid(p.avg.chin), col = 2)
library(sjPlot)
# QQ plot of random effects quantiles against standard normal quantiles
sjp.glmer(p.avg.chin, type = "re.qq")
Chum purse seine model
response = average abundance of chum salmon [per purse seine site]
This is 2nd of 4 purse models for (1)Chinook, (2)chum, (3)other migratory species and (4)resident species
Load all data - standardize variables, create subgroups from purse seine data
### Grab just chum catch - automatically removes any 0s or unidentified
### species. 2
p.2 <- purse[which(purse$Species == "Chum"), ]
### standardize variables
p.2$s.temp <- standardize(p.2$Temp.surf, centerFun = mean, scaleFun = sd)
summary(p.2$s.temp) #can see that mean is now 0, and SD is on the original scale of x (temp) -- i.e. predictor centering ## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -1.9135 -0.1349 0.1520 0.0000 0.4476 2.3575p.2$s.sal <- standardize(p.2$Sal.surf, centerFun = mean, scaleFun = sd)
p.2$s.do <- standardize(p.2$DOmg.surf, centerFun = mean, scaleFun = sd)
p.2$s.pH <- standardize(p.2$pH.surf, centerFun = mean, scaleFun = sd)
p.2$s.J.date <- standardize(p.2$J.date, centerFun = mean, scaleFun = sd)
### Create variable j2 which is Julian day squared in order to represent J
### date as quadratic relationship instead of linear
p.2$j2 <- p.2$s.J.date^2
# summarize by site-day
p.chum <- ddply(p.2, .(Year, J.date, s.J.date, j2, Habitat, Site, s.temp, s.sal,
s.do, s.pH), summarize, abundance = sum(abundance))
# plot abundance~J.date to see data distribution
plot(abundance ~ J.date, data = p.chum)
Habitat variables
Add eelgrass and sand flat habitat variables for each site
## Add habitat variables to p.chum
p.chum.hab <- p.chum
p.chum.hab$leaf_area_index <- p.hab2[match(p.chum.hab$Site, p.hab2$Site), 4]
p.chum.hab$meanturb <- p.hab2[match(p.chum.hab$Site, p.hab2$Site), 5]
## standardize
p.chum.hab$leaf_area_index <- standardize(p.chum.hab$leaf_area_index, centerFun = mean,
scaleFun = sd)
p.chum.hab$leaf_area_index <- as.numeric(p.chum.hab$leaf_area_index)
p.chum.hab$meanturb <- standardize(p.chum.hab$meanturb, centerFun = mean, scaleFun = sd)
p.chum.hab$meanturb <- as.numeric(p.chum.hab$meanturb)VIF for collinearity of habitat variables
Assess variance inflation factors
vif(lm(abundance ~ Habitat + leaf_area_index + s.J.date + j2 + meanturb + Year +
s.temp + s.sal + s.do + s.pH, data = p.chum.hab))
## Habitat leaf_area_index s.J.date j2
## 5.869989 7.362660 3.996986 2.553199
## meanturb Year s.temp s.sal
## 6.961885 1.318973 5.240322 3.951466
## s.do s.pH
## 2.674625 3.008126
# alias function identifies covariates that are multiples of each other - in
# this case we are OK.
alias(lm(abundance ~ Habitat + leaf_area_index + s.J.date + j2 + meanturb +
Year + s.temp + s.sal + s.do + s.pH, data = p.chum.hab))
## Model :
## abundance ~ Habitat + leaf_area_index + s.J.date + j2 + meanturb +
## Year + s.temp + s.sal + s.do + s.pH
## Pearson Corr with all vars
year <- p.chum.hab$Year
jday <- p.chum.hab$s.J.date
j2 <- p.chum.hab$j2
temp <- p.chum.hab$s.temp
do <- p.chum.hab$s.do
ph <- p.chum.hab$s.pH
sal <- p.chum.hab$s.sal
lai <- p.chum.hab$leaf_area_index
turb <- p.chum.hab$meanturb
hab <- p.chum.hab$Habitat
habcovar <- cbind.data.frame(hab, year, lai, turb, jday, j2, temp, do, ph, sal)
ggpairs(data = na.omit(habcovar), title = "Pearson Correlation plot habitat variables")
ggsave("/Users/Lia/Documents/Git/Fraser-salmon/Ch1.Seasonal_diversity/expl.plots/PearsonCorrAllVars_chum.pdf")Model selection
Full model
We started with a full global model with all the abiotic factors and all the habitat covariates included as well as site as a random effect. That model did not converge, so we performed manual backward selection beginning with the variables with highest VIF, and retaining the top variables from the corresponding marsh model, until the model converged. Throughout the process we viewed the VIF and model assumption plots, and used a hypothesis-based approach to select variables (as well as trying multiple iterations).
## Chum: full model with habitat variables
p.chum.1 <- glmer.nb(abundance ~ Habitat + s.J.date + j2 + Year + s.do + meanturb +
(1 | Site), data = p.chum.hab, na.action = "na.fail")
summary(p.chum.1) ## AIC 217.9## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: Negative Binomial(1.733) ( log )
## Formula: abundance ~ Habitat + s.J.date + j2 + Year + s.do + meanturb +
## (1 | Site)
## Data: p.chum.hab
##
## AIC BIC logLik deviance df.resid
## 217.9 230.5 -100.0 199.9 21
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.1241 -0.6075 -0.1823 0.4802 2.3779
##
## Random effects:
## Groups Name Variance Std.Dev.
## Site (Intercept) 1.547 1.244
## Number of obs: 30, groups: Site, 12
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.08269 0.94034 1.151 0.24957
## HabitatSand flat 0.04549 1.39273 0.033 0.97394
## s.J.date -0.98209 0.31976 -3.071 0.00213 **
## j2 -0.91609 0.23062 -3.972 0.0000712 ***
## Year 1.75525 0.63867 2.748 0.00599 **
## s.do -0.32625 0.26171 -1.247 0.21253
## meanturb -0.87781 0.65204 -1.346 0.17822
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) HbttSf s.J.dt j2 Year s.do
## HabttSndflt -0.668
## s.J.date -0.035 0.026
## j2 -0.128 -0.016 0.528
## Year -0.546 -0.033 -0.097 -0.088
## s.do -0.073 -0.186 -0.106 0.463 0.222
## meanturb 0.528 -0.780 0.207 0.134 -0.041 0.090# test VIF of full model -- several parameters are a bit high
vif.lme(p.chum.1)## HabitatSand flat s.J.date j2 Year
## 2.915721 1.914494 2.374107 1.120226
## s.do meanturb
## 1.874700 2.967015# Test fit and assumptions of full (global) model Plot residuals vs. fitted
# values
plot(fitted(p.chum.1), resid(p.chum.1), main = "Global Chum1 GLMM", xlab = "Fitted",
ylab = "Pearson residuals")
## q plot
qqnorm(x = p.chum.hab$abundance, y = resid(p.chum.1), main = "Q-Q Global Chum GLMM")
qqline(resid(p.chum.1), col = 2)
# QQ plot of random effects quantiles against standard normal quantiles --
# this is odd, all at 0???
sjp.glmer(p.chum.1, type = "re.qq")
# Check model assumptions
sjp.glmer(p.chum.1, type = "ma")
# See marginal and conditional pseudo R2 for full model - quite high.
rsquared(p.chum.1, aicc = TRUE)## Class Family Link n Marginal Conditional
## 1 glmerMod Negative Binomial log 30 0.4332422 0.9569078## Note did an iteration of global model that also converged: habitat + Jday
## + Year + DO + meanturb + Temp, AIC was higher (225.1), VIFs were higher
## and temp relationship was non linear, but marginal r squared was over 0.9Model selection - dredge
Then continued to dredge function to determine optimal model using AICc (chum has smaller sample size: 30 observations at 12 sites).
# Generate model set p.model.set.chum <- dredge(p.chum.1, extra =
# 'r.squaredGLMM') #NOTE r.squaredGLMM doesn't work with lme4. Could
# re-model with lme...
p.model.set.chum <- dredge(p.chum.1)
# Create model selection table
p.model_table.chum <- model.sel(p.model.set.chum)
options(scipen = 7)
names(p.model_table.chum) <- c("(Int)", "Habitat", "Jday^2", "Mean turbidity",
"Dissolved oxygen", "Jday", "Year", "df", "LL", "AICc", "delta", "weight")
kable(head(p.model_table.chum, n = 100), digits = 3)| (Int) | Habitat | Jday^2 | Mean turbidity | Dissolved oxygen | Jday | Year | df | LL | AICc | delta | weight | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 55 | 0.872 | NA | -0.780 | -0.908 | NA | -1.031 | 1.929 | 7 | -100.763 | 220.617 | 0.000 | 0.27914960100 |
| 51 | 0.876 | NA | -0.695 | NA | NA | -0.906 | 1.741 | 6 | -103.281 | 222.215 | 1.598 | 0.12553014162 |
| 52 | 1.658 | + | -0.767 | NA | NA | -0.975 | 1.852 | 7 | -101.577 | 222.244 | 1.627 | 0.12372240894 |
| 63 | 1.103 | NA | -0.916 | -0.861 | -0.324 | -0.983 | 1.757 | 8 | -99.963 | 222.784 | 2.167 | 0.09444822990 |
| 59 | 1.165 | NA | -0.881 | NA | -0.418 | -0.841 | 1.554 | 7 | -101.971 | 223.032 | 2.415 | 0.08343453825 |
| 56 | 1.009 | + | -0.783 | -0.802 | NA | -1.029 | 1.929 | 8 | -100.738 | 224.333 | 3.717 | 0.04352769299 |
| 60 | 1.785 | + | -0.899 | NA | -0.324 | -0.930 | 1.682 | 8 | -100.838 | 224.533 | 3.917 | 0.03938599388 |
| 27 | 2.522 | NA | -0.857 | NA | -0.533 | -0.798 | NA | 6 | -105.076 | 225.804 | 5.187 | 0.02087050651 |
| 43 | 1.130 | NA | -0.528 | NA | -0.458 | NA | 1.445 | 6 | -105.359 | 226.370 | 5.754 | 0.01571996457 |
| 23 | 2.494 | NA | -0.675 | -0.813 | NA | -0.986 | NA | 6 | -105.381 | 226.414 | 5.797 | 0.01538258219 |
| 35 | 0.803 | NA | -0.297 | NA | NA | NA | 1.672 | 5 | -107.047 | 226.594 | 5.978 | 0.01405462651 |
| 31 | 2.614 | NA | -0.883 | -0.753 | -0.470 | -0.910 | NA | 7 | -103.769 | 226.629 | 6.012 | 0.01381356356 |
| 19 | 2.367 | NA | -0.608 | NA | NA | -0.880 | NA | 5 | -107.144 | 226.788 | 6.172 | 0.01275418758 |
| 64 | 1.082 | + | -0.916 | -0.878 | -0.326 | -0.982 | 1.756 | 9 | -99.963 | 226.926 | 6.309 | 0.01190731187 |
| 33 | 0.714 | NA | NA | NA | NA | NA | 1.435 | 4 | -108.752 | 227.104 | 6.487 | 0.01089345745 |
| 39 | 0.784 | NA | -0.319 | -0.490 | NA | NA | 1.800 | 6 | -105.854 | 227.361 | 6.744 | 0.00958002895 |
| 20 | 3.140 | + | -0.667 | NA | NA | -0.944 | NA | 6 | -105.933 | 227.519 | 6.902 | 0.00885228460 |
| 28 | 3.124 | + | -0.871 | NA | -0.472 | -0.869 | NA | 7 | -104.398 | 227.887 | 7.271 | 0.00736241268 |
| 36 | 1.256 | + | -0.335 | NA | NA | NA | 1.763 | 6 | -106.129 | 227.910 | 7.293 | 0.00728029986 |
| 37 | 0.700 | NA | NA | -0.454 | NA | NA | 1.515 | 5 | -107.848 | 228.195 | 7.578 | 0.00631279883 |
| 47 | 1.077 | NA | -0.508 | -0.433 | -0.399 | NA | 1.549 | 7 | -104.627 | 228.344 | 7.728 | 0.00585855514 |
| 11 | 2.379 | NA | -0.518 | NA | -0.575 | NA | NA | 5 | -108.085 | 228.670 | 8.053 | 0.00497902565 |
| 34 | 1.056 | + | NA | NA | NA | NA | 1.470 | 5 | -108.254 | 229.009 | 8.392 | 0.00420310884 |
| 44 | 1.409 | + | -0.514 | NA | -0.399 | NA | 1.514 | 7 | -104.992 | 229.076 | 8.459 | 0.00406387588 |
| 49 | 0.691 | NA | NA | NA | NA | -0.195 | 1.402 | 5 | -108.465 | 229.429 | 8.813 | 0.00340557955 |
| 24 | 2.610 | + | -0.678 | -0.723 | NA | -0.986 | NA | 7 | -105.366 | 229.823 | 9.206 | 0.00279745513 |
| 41 | 0.738 | NA | NA | NA | -0.051 | NA | 1.396 | 5 | -108.719 | 229.937 | 9.321 | 0.00264174320 |
| 1 | 1.997 | NA | NA | NA | NA | NA | NA | 3 | -111.646 | 230.216 | 9.599 | 0.00229830952 |
| 32 | 2.510 | + | -0.883 | -0.836 | -0.476 | -0.908 | NA | 8 | -103.758 | 230.374 | 9.757 | 0.00212353075 |
| 53 | 0.673 | NA | NA | -0.529 | NA | -0.243 | 1.479 | 6 | -107.398 | 230.448 | 9.831 | 0.00204684565 |
| 40 | 0.922 | + | -0.327 | -0.387 | NA | NA | 1.803 | 7 | -105.819 | 230.730 | 10.113 | 0.00177759393 |
| 3 | 2.210 | NA | -0.217 | NA | NA | NA | NA | 4 | -110.666 | 230.931 | 10.315 | 0.00160712724 |
| 15 | 2.405 | NA | -0.506 | -0.372 | -0.545 | NA | NA | 6 | -107.656 | 230.964 | 10.347 | 0.00158144855 |
| 38 | 0.655 | + | NA | -0.490 | NA | NA | 1.516 | 6 | -107.844 | 231.340 | 10.723 | 0.00131012190 |
| 45 | 0.702 | NA | NA | -0.454 | -0.004 | NA | 1.512 | 6 | -107.847 | 231.347 | 10.730 | 0.00130567347 |
| 12 | 2.625 | + | -0.510 | NA | -0.547 | NA | NA | 6 | -107.921 | 231.493 | 10.877 | 0.00121340915 |
| 50 | 1.066 | + | NA | NA | NA | -0.208 | 1.435 | 6 | -107.926 | 231.505 | 10.888 | 0.00120636769 |
| 5 | 2.026 | NA | NA | -0.405 | NA | NA | NA | 4 | -111.037 | 231.673 | 11.057 | 0.00110890044 |
| 17 | 1.940 | NA | NA | NA | NA | -0.228 | NA | 4 | -111.225 | 232.050 | 11.433 | 0.00091868531 |
| 48 | 0.995 | + | -0.509 | -0.500 | -0.407 | NA | 1.544 | 8 | -104.616 | 232.089 | 11.473 | 0.00090071322 |
| 42 | 1.055 | + | NA | NA | 0.007 | NA | 1.477 | 6 | -108.254 | 232.160 | 11.543 | 0.00086962621 |
| 2 | 2.324 | + | NA | NA | NA | NA | NA | 4 | -111.299 | 232.199 | 11.582 | 0.00085265679 |
| 9 | 1.958 | NA | NA | NA | -0.153 | NA | NA | 4 | -111.346 | 232.291 | 11.675 | 0.00081411494 |
| 7 | 2.265 | NA | -0.228 | -0.425 | NA | NA | NA | 5 | -109.938 | 232.376 | 11.759 | 0.00078039039 |
| 57 | 0.675 | NA | NA | NA | 0.033 | -0.215 | 1.422 | 6 | -108.454 | 232.561 | 11.944 | 0.00071147652 |
| 4 | 2.638 | + | -0.240 | NA | NA | NA | NA | 5 | -110.115 | 232.731 | 12.114 | 0.00065356950 |
| 21 | 1.976 | NA | NA | -0.486 | NA | -0.269 | NA | 5 | -110.448 | 233.395 | 12.779 | 0.00046881044 |
| 61 | 0.604 | NA | NA | -0.571 | 0.129 | -0.325 | 1.569 | 7 | -107.253 | 233.597 | 12.981 | 0.00042374699 |
| 54 | 0.572 | + | NA | -0.610 | NA | -0.247 | 1.481 | 7 | -107.382 | 233.855 | 13.239 | 0.00037247648 |
| 13 | 2.002 | NA | NA | -0.390 | -0.126 | NA | NA | 5 | -110.824 | 234.148 | 13.531 | 0.00032183655 |
| 18 | 2.304 | + | NA | NA | NA | -0.245 | NA | 5 | -110.827 | 234.155 | 13.538 | 0.00032066042 |
| 16 | 2.227 | + | -0.507 | -0.514 | -0.556 | NA | NA | 7 | -107.619 | 234.329 | 13.713 | 0.00029386438 |
| 6 | 2.000 | + | NA | -0.431 | NA | NA | NA | 5 | -111.035 | 234.569 | 13.953 | 0.00026063609 |
| 58 | 1.062 | + | NA | NA | 0.129 | -0.287 | 1.523 | 7 | -107.792 | 234.675 | 14.058 | 0.00024726358 |
| 10 | 2.246 | + | NA | NA | -0.119 | NA | NA | 5 | -111.118 | 234.737 | 14.120 | 0.00023974300 |
| 46 | 0.652 | + | NA | -0.494 | -0.009 | NA | 1.508 | 7 | -107.843 | 234.777 | 14.160 | 0.00023498507 |
| 25 | 1.935 | NA | NA | NA | -0.076 | -0.176 | NA | 5 | -111.170 | 234.840 | 14.223 | 0.00022769712 |
| 8 | 2.365 | + | -0.233 | -0.350 | NA | NA | NA | 6 | -109.924 | 235.501 | 14.884 | 0.00016362014 |
| 22 | 1.881 | + | NA | -0.563 | NA | -0.272 | NA | 6 | -110.436 | 236.524 | 15.907 | 0.00009810346 |
| 29 | 1.975 | NA | NA | -0.481 | -0.014 | -0.259 | NA | 6 | -110.446 | 236.544 | 15.927 | 0.00009711378 |
| 14 | 1.861 | + | NA | -0.500 | -0.137 | NA | NA | 6 | -110.799 | 237.250 | 16.633 | 0.00006822990 |
| 26 | 2.297 | + | NA | NA | -0.017 | -0.229 | NA | 6 | -110.825 | 237.302 | 16.685 | 0.00006648867 |
| 62 | 0.579 | + | NA | -0.592 | 0.126 | -0.325 | 1.568 | 8 | -107.252 | 237.361 | 16.745 | 0.00006452894 |
| 30 | 1.865 | + | NA | -0.567 | -0.023 | -0.256 | NA | 7 | -110.431 | 239.953 | 19.337 | 0.00001765871 |
# Model averaging Version 1: use all models with delta AIC score less than 4
p.model.set.chum.4 <- get.models(p.model.set.chum, subset = delta < 4)
p.avg_model.chum.4 <- model.avg(p.model.set.chum.4)
summary(p.avg_model.chum.4)
##
## Call:
## model.avg(object = p.model.set.chum.4)
##
## Component model call:
## glmer(formula = abundance ~ <7 unique rhs>, data = p.chum.hab,
## family = negative.binomial(theta = 1.73301192961737), na.action =
## na.fail)
##
## Component models:
## df logLik AICc delta weight
## 2356 7 -100.76 220.62 0.00 0.35
## 256 6 -103.28 222.22 1.60 0.16
## 1256 7 -101.58 222.24 1.63 0.16
## 23456 8 -99.96 222.78 2.17 0.12
## 2456 7 -101.97 223.03 2.42 0.11
## 12356 8 -100.74 224.33 3.72 0.06
## 12456 8 -100.84 224.53 3.92 0.05
##
## Term codes:
## Habitat j2 meanturb s.do s.J.date Year
## 1 2 3 4 5 6
##
## Model-averaged coefficients:
## (full average)
## Estimate Std. Error Adjusted SE z value Pr(>|z|)
## (Intercept) 1.10770 0.79318 0.83059 1.334 0.182324
## j2 -0.79740 0.22508 0.23645 3.372 0.000745 ***
## meanturb -0.46845 0.53746 0.54738 0.856 0.392108
## s.J.date -0.97137 0.32964 0.34732 2.797 0.005162 **
## Year 1.81439 0.64083 0.67526 2.687 0.007210 **
## HabitatSand flat -0.34304 0.80815 0.82582 0.415 0.677851
## s.do -0.09921 0.21202 0.21707 0.457 0.647652
##
## (conditional average)
## Estimate Std. Error Adjusted SE z value Pr(>|z|)
## (Intercept) 1.1077 0.7932 0.8306 1.334 0.182324
## j2 -0.7974 0.2251 0.2364 3.372 0.000745 ***
## meanturb -0.8863 0.4197 0.4433 1.999 0.045590 *
## s.J.date -0.9714 0.3296 0.3473 2.797 0.005162 **
## Year 1.8144 0.6408 0.6753 2.687 0.007210 **
## HabitatSand flat -1.3102 1.1078 1.1565 1.133 0.257274
## s.do -0.3604 0.2630 0.2776 1.298 0.194214
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Relative variable importance:
## j2 s.J.date Year meanturb s.do Habitat
## Importance: 1.00 1.00 1.00 0.53 0.28 0.26
## N containing models: 7 7 7 3 3 3
p.chum.ci <- data.frame(confint(p.avg_model.chum.4, full = TRUE))
# Get pseudo R squared values for models up to delta < 4
p.chum.4.Rsq <- rsquared(p.model.set.chum.4, aicc = TRUE)
## write tables to .csv for easy comparison and plugging into appendix table
p.avg_model_4df.chum <- data.frame(p.avg_model.chum.4$msTable)
p.avg_model_components4.chum <- cbind(p.chum.4.Rsq, p.avg_model_4df.chum)
p.r = data.frame(Coeff = rownames(p.avg_model_4df.chum, rep(NA, length(p.avg_model_components4.chum))))
p.avg_model_components4.chum <- cbind(p.avg_model_components4.chum, p.r)
p.avg_model_components4.chum <- p.avg_model_components4.chum[, -c(7, 8)]
# write.csv(p.avg_model_components4.chum,
# '/Users/Lia/Documents/Git/Fraser-salmon/all.data/avg_model_components4_pursechum.csv')Parameter Plot
The results of model averaging including all top ranked models up to delta AICc 4
Migratory purse seine model
response = average abundance of migra salmon [per purse seine site]
This is 3rd of 4 purse models for (1)Chinook, (2)migra, (3)other migratory species and (4)resident species
Load all data - standardize variables, create subgroups from purse seine data
### Grab just migratory catch - automatically removes any 0s or unidentified
### species. 3
p.3 <- purse[which(purse$class == "migratory"), ]
p.3 <- p.3[which(p.3$Species != "Chinook"), ]
p.3 <- p.3[which(p.3$Species != "Chum"), ]
### standardize variables
p.3$s.temp <- standardize(p.3$Temp.surf, centerFun = mean, scaleFun = sd)
summary(p.3$s.temp) #can see that mean is now 0, and SD is on the original scale of x (temp) -- i.e. predictor centering ## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -3.49165 -0.69003 -0.05758 0.00000 0.38148 3.14128p.3$s.sal <- standardize(p.3$Sal.surf, centerFun = mean, scaleFun = sd)
p.3$s.do <- standardize(p.3$DOmg.surf, centerFun = mean, scaleFun = sd)
p.3$s.pH <- standardize(p.3$pH.surf, centerFun = mean, scaleFun = sd)
p.3$s.J.date <- standardize(p.3$J.date, centerFun = mean, scaleFun = sd)
### Create variable j2 which is Julian day squared in order to represent J
### date as quadratic relationship instead of linear
p.3$j2 <- p.3$s.J.date^2
# summarize by site-day
p.migra <- ddply(p.3, .(Year, J.date, s.J.date, j2, Habitat, Site, s.temp, s.sal,
s.do, s.pH), summarize, abundance = sum(abundance))
# plot abundance~J.date to see data distribution
plot(abundance ~ J.date, data = p.migra)
Habitat variables
Add eelgrass and sand flat habitat variables for each site
## Add habitat variables to p.migra
p.migra.hab <- p.migra
p.migra.hab$leaf_area_index <- p.hab2[match(p.migra.hab$Site, p.hab2$Site),
4]
p.migra.hab$meanturb <- p.hab2[match(p.migra.hab$Site, p.hab2$Site), 5]
## standardize
p.migra.hab$leaf_area_index <- standardize(p.migra.hab$leaf_area_index, centerFun = mean,
scaleFun = sd)
p.migra.hab$leaf_area_index <- as.numeric(p.migra.hab$leaf_area_index)
p.migra.hab$meanturb <- standardize(p.migra.hab$meanturb, centerFun = mean,
scaleFun = sd)
p.migra.hab$meanturb <- as.numeric(p.migra.hab$meanturb)VIF for collinearity of habitat variables
Assess variance inflation factors
vif(lm(abundance ~ Habitat + leaf_area_index + s.J.date + j2 + meanturb + Year +
s.temp + s.sal + s.do + s.pH, data = p.migra.hab))
## Habitat leaf_area_index s.J.date j2
## 5.538994 5.839467 4.459692 4.894707
## meanturb Year s.temp s.sal
## 4.645115 1.710438 3.176537 4.043365
## s.do s.pH
## 2.208030 2.635699
# alias function identifies covariates that are multiples of each other - in
# this case we are OK.
alias(lm(abundance ~ Habitat + leaf_area_index + s.J.date + j2 + meanturb +
Year + s.temp + s.sal + s.do + s.pH, data = p.migra.hab))
## Model :
## abundance ~ Habitat + leaf_area_index + s.J.date + j2 + meanturb +
## Year + s.temp + s.sal + s.do + s.pH
## Pearson Corr with all vars
year <- p.migra.hab$Year
jday <- p.migra.hab$s.J.date
j2 <- p.migra.hab$j2
temp <- p.migra.hab$s.temp
do <- p.migra.hab$s.do
ph <- p.migra.hab$s.pH
sal <- p.migra.hab$s.sal
lai <- p.migra.hab$leaf_area_index
turb <- p.migra.hab$meanturb
hab <- p.migra.hab$Habitat
habcovar <- cbind.data.frame(hab, year, lai, turb, jday, j2, temp, do, ph, sal)
ggpairs(data = na.omit(habcovar), title = "Pearson Correlation plot habitat variables")
ggsave("/Users/Lia/Documents/Git/Fraser-salmon/Ch1.Seasonal_diversity/expl.plots/PearsonCorrAllVars_migratory.pdf")Model selection
Full model
We started with a full global model with all the abiotic factors and all the habitat covariates included as well as site as a random effect. That model did not converge, so we performed manual backward selection beginning with the variables with highest VIF, and retaining the top variables from the corresponding marsh model, until the model converged. Throughout the process we viewed the VIF and model assumption plots, and used a hypothesis-based approach to select variables (as well as trying multiple iterations). In this case, the marsh model did include j^2, but looking at the purse migratory abundance distribution over time (plot above), this didn’t make sense. We removed it first and re-ran VIF and found that they all improved, supporting that this was not a good fit for this data. The marsh model also contained temp, but the distribution was quite odd with this subset and the model converged when it was removed so we left it out.
## Migratory: full model with habitat variables
p.migra.1 <- glmer.nb(abundance ~ Habitat + s.J.date + Year + s.do + s.pH +
(1 | Site), data = p.migra.hab, na.action = "na.fail")
summary(p.migra.1)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: Negative Binomial(0.6208) ( log )
## Formula: abundance ~ Habitat + s.J.date + Year + s.do + s.pH + (1 | Site)
## Data: p.migra.hab
##
## AIC BIC logLik deviance df.resid
## 722.8 742.8 -353.4 706.8 82
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -0.7721 -0.6739 -0.4504 0.1326 4.9788
##
## Random effects:
## Groups Name Variance Std.Dev.
## Site (Intercept) 0.08809 0.2968
## Number of obs: 90, groups: Site, 13
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 2.6473 0.2626 10.080 < 2e-16 ***
## HabitatSand flat 0.5550 0.3853 1.440 0.14976
## s.J.date 0.2969 0.1557 1.907 0.05654 .
## Year 0.2358 0.3525 0.669 0.50360
## s.do 0.5340 0.1701 3.139 0.00169 **
## s.pH -0.0393 0.1737 -0.226 0.82105
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) HbttSf s.J.dt Year s.do
## HabttSndflt -0.440
## s.J.date -0.157 -0.226
## Year -0.562 -0.114 0.348
## s.do 0.116 -0.322 0.452 0.070
## s.pH 0.047 0.401 -0.134 -0.335 -0.391# test VIF of full model
vif.lme(p.migra.1)## HabitatSand flat s.J.date Year s.do
## 1.261972 1.521183 1.348178 1.576720
## s.pH
## 1.541817# Test fit and assumptions of full (global) model Plot residuals vs. fitted
# values
plot(fitted(p.migra.1), resid(p.migra.1), main = "Global Migratory GLMM", xlab = "Fitted",
ylab = "Pearson residuals")
## q plot
qqnorm(x = p.migra.hab$abundance, y = resid(p.migra.1), main = "Q-Q Global Migratory GLMM")
qqline(resid(p.migra.1), col = 2)
# QQ plot of random effects quantiles against standard normal quantiles
sjp.glmer(p.migra.1, type = "re.qq")
# Check model assumptions
sjp.glmer(p.migra.1, type = "ma")
# See marginal and conditional pseudo R2 for full model
rsquared(p.migra.1, aicc = TRUE)## Class Family Link n Marginal Conditional
## 1 glmerMod Negative Binomial log 90 0.7175173 0.9011521Model selection - dredge
Then continued to dredge function to determine optimal model using AICc (migra has 90 observations at 13 sites).
# Generate model set
p.model.set.migra <- dredge(p.migra.1)
# Create model selection table
p.model_table.migra <- model.sel(p.model.set.migra)
options(scipen = 7)
names(p.model_table.migra) <- c("(Int)", "Habitat", "Dissolved oxygen", "Jday",
"pH", "Year", "df", "LL", "AICc", "delta", "weight")
kable(head(p.model_table.migra, n = 100), digits = 3)| (Int) | Habitat | Dissolved oxygen | Jday | pH | Year | df | LL | AICc | delta | weight | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| 7 | 2.967 | NA | 0.573 | 0.308 | NA | NA | 5 | -354.746 | 720.205 | 0.000 | 0.1818628771 |
| 8 | 2.749 | + | 0.526 | 0.261 | NA | NA | 6 | -353.625 | 720.262 | 0.057 | 0.1767502619 |
| 4 | 2.742 | + | 0.377 | NA | NA | NA | 5 | -355.279 | 721.271 | 1.066 | 0.1067219127 |
| 15 | 2.961 | NA | 0.600 | 0.304 | -0.090 | NA | 6 | -354.571 | 722.154 | 1.949 | 0.0686322157 |
| 24 | 2.651 | + | 0.519 | 0.292 | NA | 0.209 | 7 | -353.432 | 722.229 | 2.024 | 0.0661099024 |
| 23 | 2.889 | NA | 0.570 | 0.334 | NA | 0.165 | 6 | -354.629 | 722.269 | 2.064 | 0.0648033627 |
| 3 | 3.015 | NA | 0.392 | NA | NA | NA | 4 | -357.053 | 722.577 | 2.372 | 0.0555458697 |
| 16 | 2.749 | + | 0.526 | 0.261 | 0.001 | NA | 7 | -353.625 | 722.616 | 2.411 | 0.0544808307 |
| 12 | 2.749 | + | 0.386 | NA | -0.022 | NA | 6 | -355.268 | 723.549 | 3.343 | 0.0341780859 |
| 20 | 2.758 | + | 0.382 | NA | NA | -0.032 | 6 | -355.273 | 723.558 | 3.353 | 0.0340130156 |
| 31 | 2.835 | NA | 0.607 | 0.342 | -0.129 | 0.261 | 7 | -354.302 | 723.971 | 3.765 | 0.0276772373 |
| 11 | 3.012 | NA | 0.435 | NA | -0.122 | NA | 5 | -356.719 | 724.151 | 3.946 | 0.0252861080 |
| 32 | 2.647 | + | 0.534 | 0.297 | -0.039 | 0.236 | 8 | -353.406 | 724.590 | 4.385 | 0.0203034131 |
| 19 | 3.067 | NA | 0.406 | NA | NA | -0.111 | 5 | -356.993 | 724.701 | 4.495 | 0.0192126748 |
| 2 | 2.718 | + | NA | NA | NA | NA | 4 | -358.536 | 725.543 | 5.337 | 0.0126112635 |
| 28 | 2.757 | + | 0.387 | NA | -0.019 | -0.017 | 7 | -355.267 | 725.900 | 5.694 | 0.0105487279 |
| 27 | 3.019 | NA | 0.436 | NA | -0.120 | -0.014 | 6 | -356.718 | 726.447 | 6.242 | 0.0080220936 |
| 10 | 2.675 | + | NA | NA | 0.129 | NA | 5 | -358.185 | 727.083 | 6.878 | 0.0058372244 |
| 1 | 3.025 | NA | NA | NA | NA | NA | 3 | -360.610 | 727.499 | 7.294 | 0.0047419292 |
| 18 | 2.645 | + | NA | NA | NA | 0.148 | 5 | -358.416 | 727.546 | 7.340 | 0.0046320854 |
| 6 | 2.719 | + | NA | 0.046 | NA | NA | 5 | -358.476 | 727.666 | 7.461 | 0.0043612742 |
| 14 | 2.667 | + | NA | 0.090 | 0.166 | NA | 6 | -357.972 | 728.956 | 8.751 | 0.0022883439 |
| 26 | 2.656 | + | NA | NA | 0.118 | 0.046 | 6 | -358.175 | 729.363 | 9.158 | 0.0018672219 |
| 22 | 2.603 | + | NA | 0.093 | NA | 0.243 | 6 | -358.216 | 729.444 | 9.238 | 0.0017933045 |
| 5 | 3.020 | NA | NA | 0.061 | NA | NA | 4 | -360.517 | 729.504 | 9.298 | 0.0017402399 |
| 17 | 2.981 | NA | NA | NA | NA | 0.092 | 4 | -360.568 | 729.606 | 9.400 | 0.0016538281 |
| 9 | 3.025 | NA | NA | NA | 0.015 | NA | 4 | -360.605 | 729.680 | 9.475 | 0.0015931678 |
| 30 | 2.607 | + | NA | 0.111 | 0.139 | 0.146 | 7 | -357.892 | 731.149 | 10.944 | 0.0007644781 |
| 21 | 2.935 | NA | NA | 0.092 | NA | 0.175 | 5 | -360.387 | 731.488 | 11.283 | 0.0006451873 |
| 13 | 3.019 | NA | NA | 0.067 | 0.032 | NA | 5 | -360.496 | 731.706 | 11.501 | 0.0005785645 |
| 25 | 2.980 | NA | NA | NA | -0.002 | 0.094 | 5 | -360.567 | 731.849 | 11.644 | 0.0005386564 |
| 29 | 2.936 | NA | NA | 0.092 | 0.006 | 0.171 | 6 | -360.386 | 733.785 | 13.579 | 0.0002046418 |
# Model averaging Version 1: use all models with delta AIC score less than 4
p.model.set.migra.4 <- get.models(p.model.set.migra, subset = delta < 4)
p.avg_model.migra.4 <- model.avg(p.model.set.migra.4)
summary(p.avg_model.migra.4)
##
## Call:
## model.avg(object = p.model.set.migra.4)
##
## Component model call:
## glmer(formula = abundance ~ <12 unique rhs>, data = p.migra.hab,
## family = negative.binomial(theta = 0.620809109139615), na.action =
## na.fail)
##
## Component models:
## df logLik AICc delta weight
## 23 5 -354.75 720.21 0.00 0.20
## 123 6 -353.63 720.26 0.06 0.20
## 12 5 -355.28 721.27 1.07 0.12
## 234 6 -354.57 722.15 1.95 0.08
## 1235 7 -353.43 722.23 2.02 0.07
## 235 6 -354.63 722.27 2.06 0.07
## 2 4 -357.05 722.58 2.37 0.06
## 1234 7 -353.63 722.62 2.41 0.06
## 124 6 -355.27 723.55 3.34 0.04
## 125 6 -355.27 723.56 3.35 0.04
## 2345 7 -354.30 723.97 3.77 0.03
## 24 5 -356.72 724.15 3.95 0.03
##
## Term codes:
## Habitat s.do s.J.date s.pH Year
## 1 2 3 4 5
##
## Model-averaged coefficients:
## (full average)
## Estimate Std. Error Adjusted SE z value Pr(>|z|)
## (Intercept) 2.83840 0.24781 0.25062 11.326 < 2e-16 ***
## s.do 0.50708 0.17048 0.17243 2.941 0.00327 **
## s.J.date 0.20944 0.18325 0.18449 1.135 0.25628
## HabitatSand flat 0.33699 0.41559 0.41800 0.806 0.42013
## s.pH -0.01516 0.08198 0.08295 0.183 0.85502
## Year 0.03416 0.17355 0.17558 0.195 0.84576
##
## (conditional average)
## Estimate Std. Error Adjusted SE z value Pr(>|z|)
## (Intercept) 2.83840 0.24781 0.25062 11.326 < 2e-16 ***
## s.do 0.50708 0.17048 0.17243 2.941 0.00327 **
## s.J.date 0.29309 0.14992 0.15204 1.928 0.05388 .
## HabitatSand flat 0.63942 0.36652 0.37169 1.720 0.08538 .
## s.pH -0.06459 0.15954 0.16165 0.400 0.68945
## Year 0.15890 0.34684 0.35157 0.452 0.65128
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Relative variable importance:
## s.do s.J.date Habitat s.pH Year
## Importance: 1.00 0.71 0.53 0.23 0.21
## N containing models: 12 7 6 5 4
p.migra.ci <- data.frame(confint(p.avg_model.migra.4, full = TRUE))
# Get pseudo R squared values for models up to delta < 4 model.list.migra<-
# list(#manually list top x models from dredge)
p.migra.4.Rsq <- rsquared(p.model.set.migra.4, aicc = TRUE)
## write tables to .csv for easy comparison and plugging into appendix table
p.avg_model_4df.migra <- data.frame(p.avg_model.migra.4$msTable)
p.avg_model_components4.migra <- cbind(p.migra.4.Rsq, p.avg_model_4df.migra)
p.r = data.frame(Coeff = rownames(p.avg_model_4df.migra, rep(NA, length(p.avg_model_components4.migra))))
p.avg_model_components4.migra <- cbind(p.avg_model_components4.migra, p.r)
p.avg_model_components4.migra <- p.avg_model_components4.migra[, -c(7, 8)]
# write.csv(p.avg_model_components4.migra,
# '/Users/Lia/Documents/Git/Fraser-salmon/all.data/avg_model_components4_pursemigratory.csv')Parameter Plot
The results of model averaging including all top ranked models up to delta AICc 4
Resident purse seine model
response = average abundance of res salmon [per purse seine site]
This is 4th of 4 purse models for (1)Chinook, (2)res, (3)other migratory species and (4)resident species
Load all data - standardize variables, create subgroups from purse seine data
### Grab just resident catch - automatically removes any 0s or unidentified
### species. 4
p.4 <- purse[which(purse$class == "resident"), ]
### standardize variables
p.4$s.temp <- standardize(p.4$Temp.surf, centerFun = mean, scaleFun = sd)
summary(p.4$s.temp) #can see that mean is now 0, and SD is on the original scale of x (temp) -- i.e. predictor centering ## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -2.64937 -0.61349 0.01961 0.00000 0.55450 2.59038p.4$s.sal <- standardize(p.4$Sal.surf, centerFun = mean, scaleFun = sd)
p.4$s.do <- standardize(p.4$DOmg.surf, centerFun = mean, scaleFun = sd)
p.4$s.pH <- standardize(p.4$pH.surf, centerFun = mean, scaleFun = sd)
p.4$s.J.date <- standardize(p.4$J.date, centerFun = mean, scaleFun = sd)
### Create variable j2 which is Julian day squared in order to represent J
### date as quadratic relationship instead of linear
p.4$j2 <- p.4$s.J.date^2
# summarize by site-day
p.res <- ddply(p.4, .(Year, J.date, s.J.date, j2, Habitat, Site, s.temp, s.sal,
s.do, s.pH), summarize, abundance = sum(abundance))
# plot abundance~J.date to see data distribution
plot(abundance ~ J.date, data = p.res)
Habitat variables
Add eelgrass and sand flat habitat variables for each site
## Add habitat variables to p.res
p.res.hab <- p.res
p.res.hab$leaf_area_index <- p.hab2[match(p.res.hab$Site, p.hab2$Site), 4]
p.res.hab$meanturb <- p.hab2[match(p.res.hab$Site, p.hab2$Site), 5]
## standardize
p.res.hab$leaf_area_index <- standardize(p.res.hab$leaf_area_index, centerFun = mean,
scaleFun = sd)
p.res.hab$leaf_area_index <- as.numeric(p.res.hab$leaf_area_index)
p.res.hab$meanturb <- standardize(p.res.hab$meanturb, centerFun = mean, scaleFun = sd)
p.res.hab$meanturb <- as.numeric(p.res.hab$meanturb)VIF for collinearity of habitat variables
Assess variance inflation factors
vif(lm(abundance ~ Habitat + leaf_area_index + s.J.date + j2 + meanturb + Year +
s.temp + s.sal + s.do + s.pH, data = p.res.hab))
## Habitat leaf_area_index s.J.date j2
## 6.141154 6.308545 3.432246 4.060365
## meanturb Year s.temp s.sal
## 4.315397 1.452055 5.963063 4.107177
## s.do s.pH
## 1.935367 2.201143
# alias function identifies covariates that are multiples of each other - in
# this case we are OK.
alias(lm(abundance ~ Habitat + leaf_area_index + s.J.date + j2 + meanturb +
Year + s.temp + s.sal + s.do + s.pH, data = p.res.hab))
## Model :
## abundance ~ Habitat + leaf_area_index + s.J.date + j2 + meanturb +
## Year + s.temp + s.sal + s.do + s.pH
## Pearson Corr with all vars
year <- p.res.hab$Year
jday <- p.res.hab$s.J.date
j2 <- p.res.hab$j2
temp <- p.res.hab$s.temp
do <- p.res.hab$s.do
ph <- p.res.hab$s.pH
sal <- p.res.hab$s.sal
lai <- p.res.hab$leaf_area_index
turb <- p.res.hab$meanturb
hab <- p.res.hab$Habitat
habcovar <- cbind.data.frame(hab, year, lai, turb, jday, j2, temp, do, ph, sal)
ggpairs(data = na.omit(habcovar), title = "Pearson Correlation plot habitat variables")
ggsave("/Users/Lia/Documents/Git/Fraser-salmon/Ch1.Seasonal_diversity/expl.plots/PearsonCorrAllVars_residents.pdf")Model selection
Full model
We started with a full global model with all the abiotic factors and all the habitat covariates included as well as site as a random effect. Model did not converge - found J2 was biggest contributor (even just Year + Jday + Habitat + j2 did not converge). Once we removed that and the highly correlated habitat vars, model converged.
## Resident: full model with habitat variables
p.res.1 <- glmer.nb(abundance ~ Year + s.J.date + Habitat + s.sal + s.temp +
s.do + (1 | Site), data = p.res.hab, na.action = "na.fail")
summary(p.res.1)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: Negative Binomial(0.662) ( log )
## Formula: abundance ~ Year + s.J.date + Habitat + s.sal + s.temp + s.do +
## (1 | Site)
## Data: p.res.hab
##
## AIC BIC logLik deviance df.resid
## 2119.4 2148.5 -1050.7 2101.4 178
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -0.8062 -0.6744 -0.4913 0.1393 8.9313
##
## Random effects:
## Groups Name Variance Std.Dev.
## Site (Intercept) 0.1537 0.392
## Number of obs: 187, groups: Site, 13
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 5.9729 0.2325 25.692 < 2e-16 ***
## Year -0.1753 0.2244 -0.781 0.435
## s.J.date 0.2195 0.1254 1.751 0.080 .
## HabitatSand flat -2.2840 0.3808 -5.998 0.00000000199 ***
## s.sal -0.1503 0.1695 -0.887 0.375
## s.temp 0.7147 0.1395 5.125 0.00000029745 ***
## s.do 0.1235 0.1160 1.065 0.287
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) Year s.J.dt HbttSf s.sal s.temp
## Year -0.375
## s.J.date -0.007 0.152
## HabttSndflt -0.664 0.024 -0.046
## s.sal -0.386 0.235 -0.126 0.632
## s.temp -0.150 0.332 -0.458 0.193 0.429
## s.do -0.047 0.053 0.374 0.122 0.237 -0.220# test VIF of full model
vif.lme(p.res.1)## Year s.J.date HabitatSand flat s.sal
## 1.332541 1.601720 1.736491 2.321479
## s.temp s.do
## 1.904814 1.361182# Test fit and assumptions of full (global) model Plot residuals vs. fitted
# values
plot(fitted(p.res.1), resid(p.res.1), main = "Global Resident GLMM", xlab = "Fitted",
ylab = "Pearson residuals")
## q plot
qqnorm(x = p.res.hab$abundance, y = resid(p.res.1), main = "Q-Q Global Resident GLMM")
qqline(resid(p.res.1), col = 2)
# QQ plot of random effects quantiles against standard normal quantiles
sjp.glmer(p.res.1, type = "re.qq")
# Check model assumptions
sjp.glmer(p.res.1, type = "ma")
# See marginal and conditional pseudo R2 for full model -- very high
rsquared(p.res.1, aicc = TRUE)## Class Family Link n Marginal Conditional
## 1 glmerMod Negative Binomial log 187 0.9187988 0.9962039Model selection - dredge
Then continued to dredge function to determine optimal model using AICc (probably unnecessary but have done for all others: res has 187 observations at 13 sites).
# Generate model set
p.model.set.res <- dredge(p.res.1)
# Create model selection table
p.model_table.res <- model.sel(p.model.set.res)
options(scipen = 7)
names(p.model_table.res) <- c("(Int)", "Habitat", "Dissolved oxygen", "Jday",
"Salinity", "Temp", "Year", "df", "LL", "AICc", "delta", "weight")
kable(head(p.model_table.res, n = 100), digits = 3)| (Int) | Habitat | Dissolved oxygen | Jday | Salinity | Temp | Year | df | LL | AICc | delta | weight | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 18 | 5.863 | + | NA | NA | NA | 0.908 | NA | 5 | -1053.126 | 2116.584 | 0.000 | 1.531220e-01 |
| 22 | 5.848 | + | NA | 0.148 | NA | 0.848 | NA | 6 | -1052.116 | 2116.699 | 0.115 | 1.445893e-01 |
| 24 | 5.857 | + | 0.148 | 0.218 | NA | 0.787 | NA | 7 | -1051.291 | 2117.207 | 0.622 | 1.121679e-01 |
| 30 | 5.918 | + | NA | 0.181 | -0.161 | 0.788 | NA | 7 | -1051.618 | 2117.861 | 1.277 | 8.087970e-02 |
| 50 | 5.917 | + | NA | NA | NA | 0.870 | -0.188 | 6 | -1052.740 | 2117.946 | 1.362 | 7.748999e-02 |
| 26 | 5.901 | + | NA | NA | -0.078 | 0.885 | NA | 6 | -1052.996 | 2118.458 | 1.874 | 5.999424e-02 |
| 20 | 5.869 | + | 0.045 | NA | NA | 0.899 | NA | 6 | -1053.030 | 2118.526 | 1.942 | 5.798267e-02 |
| 54 | 5.884 | + | NA | 0.136 | NA | 0.827 | -0.123 | 7 | -1051.957 | 2118.540 | 1.956 | 5.758843e-02 |
| 32 | 5.908 | + | 0.128 | 0.233 | -0.120 | 0.752 | NA | 8 | -1051.026 | 2118.860 | 2.276 | 4.907511e-02 |
| 56 | 5.895 | + | 0.148 | 0.206 | NA | 0.767 | -0.127 | 8 | -1051.119 | 2119.047 | 2.463 | 4.469397e-02 |
| 62 | 5.987 | + | NA | 0.168 | -0.191 | 0.746 | -0.185 | 8 | -1051.284 | 2119.378 | 2.793 | 3.788309e-02 |
| 58 | 5.997 | + | NA | NA | -0.131 | 0.820 | -0.247 | 7 | -1052.401 | 2119.428 | 2.843 | 3.694746e-02 |
| 52 | 5.929 | + | 0.060 | NA | NA | 0.855 | -0.206 | 7 | -1052.574 | 2119.773 | 3.189 | 3.108460e-02 |
| 64 | 5.973 | + | 0.124 | 0.220 | -0.150 | 0.715 | -0.175 | 9 | -1050.722 | 2120.462 | 3.877 | 2.203178e-02 |
| 28 | 5.897 | + | 0.031 | NA | -0.062 | 0.884 | NA | 7 | -1052.956 | 2120.537 | 3.953 | 2.121642e-02 |
| 60 | 5.994 | + | 0.037 | NA | -0.114 | 0.818 | -0.251 | 8 | -1052.341 | 2121.492 | 4.908 | 1.316368e-02 |
| 17 | 4.848 | NA | NA | NA | NA | 0.909 | NA | 4 | -1063.257 | 2134.734 | 18.149 | 1.753690e-05 |
| 21 | 4.850 | NA | NA | 0.138 | NA | 0.855 | NA | 5 | -1062.396 | 2135.123 | 18.539 | 1.443497e-05 |
| 49 | 4.928 | NA | NA | NA | NA | 0.876 | -0.191 | 5 | -1062.888 | 2136.108 | 19.524 | 8.819499e-06 |
| 23 | 4.845 | NA | 0.103 | 0.188 | NA | 0.809 | NA | 6 | -1062.012 | 2136.490 | 19.906 | 7.287296e-06 |
| 25 | 4.857 | NA | NA | NA | 0.035 | 0.919 | NA | 5 | -1063.235 | 2136.801 | 20.217 | 6.237934e-06 |
| 19 | 4.847 | NA | 0.011 | NA | NA | 0.907 | NA | 5 | -1063.251 | 2136.833 | 20.249 | 6.138807e-06 |
| 53 | 4.905 | NA | NA | 0.126 | NA | 0.836 | -0.132 | 6 | -1062.226 | 2136.919 | 20.335 | 5.880206e-06 |
| 29 | 4.838 | NA | NA | 0.148 | -0.046 | 0.838 | NA | 6 | -1062.361 | 2137.188 | 20.604 | 5.139473e-06 |
| 51 | 4.929 | NA | 0.023 | NA | NA | 0.869 | -0.197 | 6 | -1062.864 | 2138.194 | 21.610 | 3.108799e-06 |
| 57 | 4.927 | NA | NA | NA | -0.012 | 0.871 | -0.196 | 6 | -1062.886 | 2138.238 | 21.654 | 3.040038e-06 |
| 55 | 4.900 | NA | 0.102 | 0.176 | NA | 0.791 | -0.132 | 7 | -1061.841 | 2138.309 | 21.724 | 2.935337e-06 |
| 31 | 4.843 | NA | 0.101 | 0.190 | -0.009 | 0.807 | NA | 7 | -1062.010 | 2138.646 | 22.062 | 2.479189e-06 |
| 27 | 4.859 | NA | 0.023 | NA | 0.049 | 0.918 | NA | 6 | -1063.214 | 2138.894 | 22.309 | 2.190767e-06 |
| 61 | 4.896 | NA | NA | 0.140 | -0.077 | 0.804 | -0.157 | 7 | -1062.137 | 2138.899 | 22.315 | 2.184566e-06 |
| 59 | 4.930 | NA | 0.023 | NA | 0.001 | 0.870 | -0.197 | 7 | -1062.864 | 2140.353 | 23.768 | 1.056271e-06 |
| 63 | 4.896 | NA | 0.096 | 0.180 | -0.040 | 0.778 | -0.145 | 8 | -1061.819 | 2140.447 | 23.863 | 1.007740e-06 |
| 48 | 6.232 | + | 0.216 | 0.608 | -0.555 | NA | -0.470 | 8 | -1064.164 | 2145.136 | 28.552 | 9.661352e-08 |
| 46 | 6.291 | + | NA | 0.509 | -0.665 | NA | -0.525 | 7 | -1065.749 | 2146.124 | 29.540 | 5.895956e-08 |
| 16 | 6.058 | + | 0.266 | 0.691 | -0.511 | NA | NA | 7 | -1066.732 | 2148.090 | 31.506 | 2.206155e-08 |
| 14 | 6.111 | + | NA | 0.587 | -0.654 | NA | NA | 6 | -1068.989 | 2150.445 | 33.861 | 6.796668e-09 |
| 40 | 5.961 | + | 0.382 | 0.679 | NA | NA | -0.357 | 7 | -1069.673 | 2153.972 | 37.388 | 1.165206e-09 |
| 8 | 5.833 | + | 0.418 | 0.726 | NA | NA | NA | 6 | -1071.259 | 2154.984 | 38.400 | 7.025299e-10 |
| 42 | 6.396 | + | NA | NA | -0.632 | NA | -0.746 | 6 | -1074.538 | 2161.543 | 44.959 | 2.644528e-11 |
| 38 | 5.972 | + | NA | 0.499 | NA | NA | -0.444 | 6 | -1074.977 | 2162.421 | 45.837 | 1.704984e-11 |
| 44 | 6.405 | + | -0.048 | NA | -0.656 | NA | -0.749 | 7 | -1074.450 | 2163.526 | 46.941 | 9.814420e-12 |
| 6 | 5.804 | + | NA | 0.544 | NA | NA | NA | 5 | -1077.488 | 2165.307 | 48.723 | 4.027507e-12 |
| 47 | 4.965 | NA | 0.212 | 0.620 | -0.499 | NA | -0.467 | 7 | -1076.309 | 2167.243 | 50.659 | 1.529650e-12 |
| 45 | 4.962 | NA | NA | 0.525 | -0.619 | NA | -0.526 | 6 | -1077.765 | 2167.996 | 51.412 | 1.049658e-12 |
| 15 | 4.787 | NA | 0.267 | 0.699 | -0.428 | NA | NA | 6 | -1078.618 | 2169.704 | 53.119 | 4.469977e-13 |
| 13 | 4.752 | NA | NA | 0.599 | -0.585 | NA | NA | 5 | -1080.788 | 2171.907 | 55.322 | 1.485640e-13 |
| 39 | 5.037 | NA | 0.366 | 0.680 | NA | NA | -0.331 | 6 | -1080.071 | 2172.608 | 56.024 | 1.046052e-13 |
| 10 | 6.097 | + | NA | NA | -0.533 | NA | NA | 5 | -1081.143 | 2172.618 | 56.033 | 1.041219e-13 |
| 7 | 4.896 | NA | 0.395 | 0.723 | NA | NA | NA | 5 | -1081.323 | 2172.978 | 56.394 | 8.695593e-14 |
| 12 | 6.103 | + | -0.029 | NA | -0.551 | NA | NA | 6 | -1081.115 | 2174.697 | 58.113 | 3.680686e-14 |
| 34 | 6.030 | + | NA | NA | NA | NA | -0.567 | 5 | -1082.221 | 2174.773 | 58.189 | 3.544188e-14 |
| 36 | 6.033 | + | 0.114 | NA | NA | NA | -0.561 | 6 | -1081.721 | 2175.908 | 59.324 | 2.009358e-14 |
| 37 | 5.077 | NA | NA | 0.523 | NA | NA | -0.410 | 5 | -1084.888 | 2180.108 | 63.523 | 2.460938e-15 |
| 2 | 5.824 | + | NA | NA | NA | NA | NA | 4 | -1086.409 | 2181.037 | 64.453 | 1.546059e-15 |
| 4 | 5.834 | + | 0.129 | NA | NA | NA | NA | 5 | -1085.779 | 2181.889 | 65.305 | 1.009782e-15 |
| 5 | 4.903 | NA | NA | 0.566 | NA | NA | NA | 4 | -1086.866 | 2181.952 | 65.368 | 9.784846e-16 |
| 41 | 5.058 | NA | NA | NA | -0.579 | NA | -0.745 | 5 | -1086.620 | 2183.571 | 66.986 | 4.356292e-16 |
| 43 | 5.052 | NA | -0.055 | NA | -0.611 | NA | -0.749 | 6 | -1086.511 | 2185.488 | 68.904 | 1.670318e-16 |
| 33 | 5.122 | NA | NA | NA | NA | NA | -0.559 | 4 | -1092.449 | 2193.118 | 76.534 | 3.680100e-18 |
| 9 | 4.771 | NA | NA | NA | -0.455 | NA | NA | 4 | -1092.671 | 2193.562 | 76.978 | 2.948248e-18 |
| 35 | 5.119 | NA | 0.111 | NA | NA | NA | -0.554 | 5 | -1091.977 | 2194.285 | 77.701 | 2.053822e-18 |
| 11 | 4.768 | NA | -0.022 | NA | -0.471 | NA | NA | 5 | -1092.655 | 2195.641 | 79.057 | 1.042326e-18 |
| 1 | 4.884 | NA | NA | NA | NA | NA | NA | 3 | -1096.177 | 2198.486 | 81.902 | 2.513613e-19 |
| 3 | 4.883 | NA | 0.121 | NA | NA | NA | NA | 4 | -1095.621 | 2199.462 | 82.878 | 1.542993e-19 |
# Model averaging Version 1: use all models with delta AIC score less than 4
p.model.set.res.4 <- get.models(p.model.set.res, subset = delta < 4)
p.avg_model.res.4 <- model.avg(p.model.set.res.4)
summary(p.avg_model.res.4)
##
## Call:
## model.avg(object = p.model.set.res.4)
##
## Component model call:
## glmer(formula = abundance ~ <15 unique rhs>, data = p.res.hab,
## family = negative.binomial(theta = 0.662012821216055), na.action =
## na.fail)
##
## Component models:
## df logLik AICc delta weight
## 15 5 -1053.13 2116.58 0.00 0.16
## 135 6 -1052.12 2116.70 0.11 0.15
## 1235 7 -1051.29 2117.21 0.62 0.11
## 1345 7 -1051.62 2117.86 1.28 0.08
## 156 6 -1052.74 2117.95 1.36 0.08
## 145 6 -1053.00 2118.46 1.87 0.06
## 125 6 -1053.03 2118.53 1.94 0.06
## 1356 7 -1051.96 2118.54 1.96 0.06
## 12345 8 -1051.03 2118.86 2.28 0.05
## 12356 8 -1051.12 2119.05 2.46 0.05
## 13456 8 -1051.28 2119.38 2.79 0.04
## 1456 7 -1052.40 2119.43 2.84 0.04
## 1256 7 -1052.57 2119.77 3.19 0.03
## 123456 9 -1050.72 2120.46 3.88 0.02
## 1245 7 -1052.96 2120.54 3.95 0.02
##
## Term codes:
## Habitat s.do s.J.date s.sal s.temp Year
## 1 2 3 4 5 6
##
## Model-averaged coefficients:
## (full average)
## Estimate Std. Error Adjusted SE z value Pr(>|z|)
## (Intercept) 5.89148 0.21936 0.22079 26.684 <2e-16 ***
## HabitatSand flat -2.15828 0.34126 0.34339 6.285 <2e-16 ***
## s.temp 0.83582 0.13190 0.13263 6.302 <2e-16 ***
## s.J.date 0.10158 0.12777 0.12817 0.793 0.428
## s.do 0.03780 0.08770 0.08803 0.429 0.668
## s.sal -0.04085 0.11108 0.11158 0.366 0.714
## Year -0.05453 0.14802 0.14871 0.367 0.714
##
## (conditional average)
## Estimate Std. Error Adjusted SE z value Pr(>|z|)
## (Intercept) 5.8915 0.2194 0.2208 26.684 <2e-16 ***
## HabitatSand flat -2.1583 0.3413 0.3434 6.285 <2e-16 ***
## s.temp 0.8358 0.1319 0.1326 6.302 <2e-16 ***
## s.J.date 0.1826 0.1206 0.1214 1.504 0.132
## s.do 0.1103 0.1202 0.1209 0.912 0.362
## s.sal -0.1308 0.1666 0.1676 0.781 0.435
## Year -0.1749 0.2219 0.2233 0.783 0.434
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Relative variable importance:
## Habitat s.temp s.J.date s.do s.sal Year
## Importance: 1.00 1.00 0.56 0.34 0.31 0.31
## N containing models: 15 15 8 7 7 7
p.res.ci <- data.frame(confint(p.avg_model.res.4, full = TRUE))
# Get pseudo R squared values for models up to delta < 4 model.list.res<-
# list(#manually list top x models from dredge)
p.res.4.Rsq <- rsquared(p.model.set.res.4, aicc = TRUE)
## write tables to .csv for easy comparison and plugging into appendix table
p.avg_model_4df.res <- data.frame(p.avg_model.res.4$msTable)
p.avg_model_components4.res <- cbind(p.res.4.Rsq, p.avg_model_4df.res)
p.r = data.frame(Coeff = rownames(p.avg_model_4df.res, rep(NA, length(p.avg_model_components4.res))))
p.avg_model_components4.res <- cbind(p.avg_model_components4.res, p.r)
p.avg_model_components4.res <- p.avg_model_components4.res[, -c(7, 8)]
# write.csv(p.avg_model_components4.res,
# '/Users/Lia/Documents/Git/Fraser-salmon/all.data/avg_model_components4_purseresidents.csv')Parameter Plot
The results of model averaging including all top ranked models up to delta AICc 4
## List of 11
## $ axis.title.x : list()
## ..- attr(*, "class")= chr [1:2] "element_blank" "element"
## $ axis.title.y :List of 11
## ..$ family : NULL
## ..$ face : NULL
## ..$ colour : NULL
## ..$ size : num 14
## ..$ hjust : NULL
## ..$ vjust : NULL
## ..$ angle : num 90
## ..$ lineheight : NULL
## ..$ margin : NULL
## ..$ debug : NULL
## ..$ inherit.blank: logi FALSE
## ..- attr(*, "class")= chr [1:2] "element_text" "element"
## $ axis.text.x : list()
## ..- attr(*, "class")= chr [1:2] "element_blank" "element"
## $ axis.text.y :List of 11
## ..$ family : NULL
## ..$ face : NULL
## ..$ colour : NULL
## ..$ size : num 12
## ..$ hjust : NULL
## ..$ vjust : NULL
## ..$ angle : NULL
## ..$ lineheight : NULL
## ..$ margin : NULL
## ..$ debug : NULL
## ..$ inherit.blank: logi FALSE
## ..- attr(*, "class")= chr [1:2] "element_text" "element"
## $ axis.ticks.x : list()
## ..- attr(*, "class")= chr [1:2] "element_blank" "element"
## $ axis.line.x :List of 6
## ..$ colour : chr "black"
## ..$ size : num 0.5
## ..$ linetype : NULL
## ..$ lineend : NULL
## ..$ arrow : logi FALSE
## ..$ inherit.blank: logi FALSE
## ..- attr(*, "class")= chr [1:2] "element_line" "element"
## $ axis.line.y :List of 6
## ..$ colour : chr "black"
## ..$ size : num 0.5
## ..$ linetype : NULL
## ..$ lineend : NULL
## ..$ arrow : logi FALSE
## ..$ inherit.blank: logi FALSE
## ..- attr(*, "class")= chr [1:2] "element_line" "element"
## $ panel.border : list()
## ..- attr(*, "class")= chr [1:2] "element_blank" "element"
## $ panel.grid.major: list()
## ..- attr(*, "class")= chr [1:2] "element_blank" "element"
## $ panel.grid.minor: list()
## ..- attr(*, "class")= chr [1:2] "element_blank" "element"
## $ plot.margin :Classes 'margin', 'unit' atomic [1:4] 2 2 2 2
## .. ..- attr(*, "valid.unit")= int 8
## .. ..- attr(*, "unit")= chr "pt"
## - attr(*, "class")= chr [1:2] "theme" "gg"
## - attr(*, "complete")= logi FALSE
## - attr(*, "validate")= logi TRUE